Merge pull request #64 from bska/ert

Bring ERT branch up-to-date with recent changes in master

Fix merge conflicts in the process.

.gitignore

@@ -1,6 +1,7 @@
 *~
 .\#*
 *.o
+*.mod
 *.lo
 *.la
 .libs
@@ -34,6 +35,7 @@ tutorials/tutorial[1-4]
 
 # Ignoring executables
 *_test
+examples/compute_tof
 examples/scaneclipsedeck
 examples/spu_2p
 examples/reorder-qfs
@@ -42,8 +44,10 @@ examples/sim_2p_incomp_reorder
 examples/sim_2p_comp_reorder
 examples/sim_wateroil
 examples/wells_example
+tests/test_agmg
 tests/test_cfs_tpfa
 tests/test_jacsys
+tests/test_read_grid
 tests/test_readvector
 tests/test_sf2p
 tests/bo_fluid_p_and_z_deps
@@ -53,4 +57,5 @@ tests/test_column_extract
 tests/test_lapack
 tests/test_read_vag
 tests/test_readpolymer
+tests/test_wells
 tests/test_writeVtkData

Makefile.am

@@ -13,7 +13,8 @@ lib_LTLIBRARIES = lib/libopmcore.la
 
 AM_CPPFLAGS =					\
 $(ERT_CPPFLAGS)                                 \
-$(OPM_BOOST_CPPFLAGS)
+$(OPM_BOOST_CPPFLAGS)                           \
+${SUPERLU_CPPFLAGS}
 
 # ----------------------------------------------------------------------
 # Link-time flags needed both to successfully link the library and to
@@ -23,6 +24,7 @@ lib_libopmcore_la_LDFLAGS =			\
 -R $(OPM_BOOST_LIBDIR)				\
 $(OPM_BOOST_LDFLAGS)                            \
 $(ERT_LDFLAGS)
+${SUPERLU_LDFLAGS}
 
 lib_libopmcore_la_LIBADD = 			\
 $(BOOST_FILESYSTEM_LIB)				\
@@ -30,7 +32,7 @@ $(BOOST_SYSTEM_LIB)				\
 $(BOOST_DATE_TIME_LIB)				\
 $(BOOST_UNIT_TEST_FRAMEWORK_LIB)		\
 $(ERT_LIBS)                                     \
-$(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS)
+$(LAPACK_LIBS) ${SUPERLU_LIBS} $(BLAS_LIBS) $(LIBS)
 
 # ----------------------------------------------------------------------
 # Library constituents.  SOURCES followed by HEADERS.
@@ -108,6 +110,8 @@ opm/core/simulator/SimulatorReport.cpp                  \
 opm/core/simulator/SimulatorTimer.cpp			\
 opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp	\
 opm/core/transport/reorder/TransportModelInterface.cpp	\
+opm/core/transport/reorder/TransportModelTracerTof.cpp	\
+opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp	\
 opm/core/transport/reorder/TransportModelTwophase.cpp	\
 opm/core/transport/reorder/nlsolvers.c			\
 opm/core/transport/reorder/reordersequence.cpp		\
@@ -230,6 +234,8 @@ opm/core/transport/SimpleFluid2pWrapper.hpp             \
 opm/core/transport/SinglePointUpwindTwoPhase.hpp	\
 opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp	\
 opm/core/transport/reorder/TransportModelInterface.hpp	\
+opm/core/transport/reorder/TransportModelTracerTof.hpp	\
+opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp	\
 opm/core/transport/reorder/TransportModelTwophase.hpp	\
 opm/core/transport/reorder/nlsolvers.h			\
 opm/core/transport/reorder/reordersequence.h		\
@@ -289,6 +295,7 @@ opm/core/linalg/call_umfpack.h                  \
 opm/core/linalg/LinearSolverUmfpack.hpp
 endif
 
+
 if HAVE_ERT
 lib_libopmcore_la_SOURCES  +=			\
 opm/core/utility/writeECLData.cpp

README

@@ -56,10 +56,10 @@ sudo add-apt-repository -y ppa:opm/ppa
 sudo apt-get update
 
 # parts of DUNE needed
-sudo apt-get install libdune-common-dev libdune-istl-dev
+sudo apt-get install libdune-common-dev libdune-istl-dev libdune-grid-dev
 
 # libraries necessary for OPM
-sudo apt-get install -y libxml0-dev
+sudo apt-get install -y libxml2-dev
 
 
 DEPENDENCIES FOR SUSE BASED DISTRIBUTIONS
@@ -88,10 +88,25 @@ If you want to contribute, fork OPM/opm-core on github.
 BUILDING
 --------
 
+There are two ways to build the opm-core library:
+
+1. As a stand-alone library.
  cd opm-core
  autoreconf -i
  ./configure
  make
+If you want to install the library:
+ make install
+or (if installing to /usr/local or similar)
+ sudo make install
+
+2. As a dune module.
+ - Put the opm-core directory in the same directory
+   as the other dune modules to be built (e.g. dune-commmon,
+   dune-grid).
+ - Run dunecontrol as normal. For more information on
+   the dune build system, see
+   http://www.dune-project.org/doc/installation-notes.html
 
 
 DOCUMENTATION

examples/Makefile.am

@@ -3,19 +3,34 @@ ERT_INCLUDE_PATH = $(ERT_ROOT)/include
 
 AM_CPPFLAGS =					\
 -I$(top_srcdir)					\
-$(OPM_BOOST_CPPFLAGS)				\
+-I$(ERT_INCLUDE_PATH)				\
--I$(ERT_INCLUDE_PATH)
+$(OPM_BOOST_CPPFLAGS)
 
 # All targets link to the library
 LDADD =						\
 $(top_builddir)/lib/libopmcore.la
 
+# Convenience definition for targets that use Boost.Filesystem directly.
+# While libopmcore depends on (and references) Boost.Filesystem (through
+# the $(BOOST_FILESYSTEM_LIB) macro) this indirect dependency is not
+# sufficient to satisfy the requirements of targets that use the indirect
+# libraries directly.
+#
+# Additional details at
+#   https://fedoraproject.org/wiki/UnderstandingDSOLinkChange
+#
+LINK_BOOST_FILESYSTEM =				\
+$(OPM_BOOST_LDFLAGS)				\
+$(BOOST_FILESYSTEM_LIB)				\
+$(BOOST_SYSTEM_LIB)
+
 # ----------------------------------------------------------------------
 # Declare products (i.e., the example programs).
 #
 # Please keep the list sorted.
 
 noinst_PROGRAMS =				\
+compute_tof					\
 refine_wells					\
 scaneclipsedeck					\
 sim_2p_comp_reorder				\
@@ -30,10 +45,19 @@ wells_example
 #
 # Please maintain sort order from "noinst_PROGRAMS".
 
+compute_tof_SOURCES           = compute_tof.cpp
+compute_tof_LDADD             = $(LDADD) $(LINK_BOOST_FILESYSTEM)
 refine_wells_SOURCES          = refine_wells.cpp
+
 sim_2p_comp_reorder_SOURCES   = sim_2p_comp_reorder.cpp
+sim_2p_comp_reorder_LDADD     = $(LDADD) $(LINK_BOOST_FILESYSTEM)
+
 sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp
+sim_2p_incomp_reorder_LDADD   = $(LDADD) $(LINK_BOOST_FILESYSTEM)
+
 sim_wateroil_SOURCES          = sim_wateroil.cpp
+sim_wateroil_LDADD            = $(LDADD) $(LINK_BOOST_FILESYSTEM)
+
 wells_example_SOURCES         = wells_example.cpp
 
 # ----------------------------------------------------------------------
@@ -45,5 +69,6 @@ noinst_PROGRAMS += spu_2p
 spu_2p_SOURCES   = spu_2p.cpp
 spu_2p_LDADD     =				\
 $(LDADD)					\
+$(LINK_BOOST_FILESYSTEM)			\
 $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS)
 endif

examples/compute_tof.cpp

@@ -0,0 +1,254 @@
+/*
+  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/utility/StopWatch.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/linalg/LinearSolverFactory.hpp>
+
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/pressure/IncompTpfa.hpp>
+#include <opm/core/transport/reorder/TransportModelTracerTof.hpp>
+#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.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 incompressible tof computations     ===============\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<IncompPropertiesInterface> props;
+    boost::scoped_ptr<Opm::WellsManager> wells;
+    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");
+        deck.reset(new EclipseGridParser(deck_filename));
+        // Grid init
+        grid.reset(new GridManager(*deck));
+        // Rock and fluid init
+        props.reset(new IncompPropertiesFromDeck(*deck, *grid->c_grid()));
+        // Wells init.
+        wells.reset(new Opm::WellsManager(*deck, *grid->c_grid(), props->permeability()));
+        // 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);
+        }
+    } 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 IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
+        // Wells init.
+        wells.reset(new Opm::WellsManager());
+        // 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;
+        }
+    }
+    const double *grav = use_gravity ? &gravity[0] : 0;
+
+    // Initialising src
+    std::vector<double> porevol;
+    computePorevolume(*grid->c_grid(), props->porosity(), porevol);
+    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 {
+        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);
+
+    // Pressure solver.
+    Opm::IncompTpfa psolver(*grid->c_grid(), *props, 0, linsolver,
+                            0.0, 0.0, 0,
+                            grav, wells->c_wells(), src, bcs.c_bcs());
+
+    // Choice of tof solver.
+    bool use_dg = param.getDefault("use_dg", false);
+    int dg_degree = -1;
+    if (use_dg) {
+        dg_degree = param.getDefault("dg_degree", 0);
+    }
+
+    // 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");
+    }
+
+    // Init wells.
+    Opm::WellState well_state;
+    well_state.init(wells->c_wells(), state);
+
+    // Main solvers.
+    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 solvers     ===============" << std::endl;
+
+    // Solve pressure.
+    pressure_timer.start();
+    psolver.solve(1.0, state, well_state);
+    pressure_timer.stop();
+    double pt = pressure_timer.secsSinceStart();
+    std::cout << "Pressure solver took:  " << pt << " seconds." << std::endl;
+    ptime += pt;
+
+    // Process transport sources (to include bdy terms and well flows).
+    std::vector<double> transport_src;
+    Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
+                                wells->c_wells(), well_state.perfRates(), transport_src);
+
+    // Solve time-of-flight.
+    std::vector<double> tof;
+    if (use_dg) {
+        Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid());
+        transport_timer.start();
+        tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof);
+        transport_timer.stop();
+    } else {
+        Opm::TransportModelTracerTof tofsolver(*grid->c_grid());
+        transport_timer.start();
+        tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof);
+        transport_timer.stop();
+    }
+    double tt = transport_timer.secsSinceStart();
+    std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
+    ttime += tt;
+    total_timer.stop();
+
+    // Output.
+    if (output) {
+        std::string tof_filename = output_dir + "/tof.txt";
+        std::ofstream tof_stream(tof_filename.c_str());
+        std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tof_stream, "\n"));
+    }
+
+    std::cout << "\n\n================    End of simulation     ===============\n"
+              << "Total time taken: " << total_timer.secsSinceStart()
+              << "\n  Pressure time:  " << ptime
+              << "\n  Transport time: " << ttime << std::endl;
+}

m4/opm_core.m4

@@ -17,8 +17,8 @@ AX_BOOST_SYSTEM
 AX_BOOST_DATE_TIME
 AX_BOOST_FILESYSTEM
 AX_BOOST_UNIT_TEST_FRAMEWORK
-
 AX_DUNE_ISTL
+OPM_PATH_SUPERLU
 OPM_AGMG
 
 # Checks for header files.

m4/opm_superlu.m4

@@ -0,0 +1,334 @@
+## -*- autoconf -*-
+# $Id: superlu.m4 5908 2010-02-23 16:46:05Z joe $
+# searches for SuperLU headers and libs
+
+# _slu_lib_path(SUPERLU_ROOT, HEADER)
+#
+# Try to find the subpath unter SUPERLU_ROOT containing HEADER.  Try
+# SUPERLU_ROOT/"include/superlu", SUPERLU_ROOT/"include", and
+# SUPERLU_ROOT/"SRC", in that order.  Set the subpath for the library to
+# "lib".  If HEADER was found in SUPERLU_ROOT/"SRC", check whether
+# SUPERLU_ROOT/"lib" is a directory, and set the subpath for the library to
+# the empty string "" if it isn't.
+#
+# Shell variables:
+#   my_include_path
+#      The subpath HEADER was found in: "include/superlu", "include", or
+#      "SRC".  Contents is only meaningful for my_slu_found=yes.
+#   my_lib_path
+#      The subpath for the library: "lib" or "". Contents is only meaningful
+#      for my_slu_found=yes.
+#   my_slu_found
+#      Whether HEADER was found at all.  Either "yes" or "no".
+AC_DEFUN([_slu_lib_path],
+    [
+        my_include_path=include/superlu
+        my_lib_path=lib
+        my_slu_found=yes
+        if test ! -f "$1/$my_include_path/$2" ; then
+            #Try to find headers under superlu
+            my_include_path=include
+            if test ! -f "$1/$my_include_path/$2" ; then
+                my_include_path=SRC
+                if test ! -f "$1/$my_include_path/$2"; then
+                    my_slu_found=no
+                else
+                    if ! test -d "$1/$my_lib_path"; then
+                        my_lib_path=""
+                    fi
+                fi
+            fi
+        fi
+    ]
+)
+
+# _slu_search_versions(SUPERLU_ROOT)
+#
+# Search for either "slu_ddefs.h" or "dsp_defs.h" using _slu_lib_path().
+#
+# Shell variables:
+#   my_slu_header
+#      The name of the header that was found: first of "slu_ddefs.h" or
+#      "dsp_defs.h".  Contents is only meaningful for my_slu_found=yes.
+#   my_include_path
+#      The subpath the header was found in: "include/superlu", "include", or
+#      "SRC".  Contents is only meaningful for my_slu_found=yes.
+#   my_lib_path
+#      The subpath for the library: "lib" or "". Contents is only meaningful
+#      for my_slu_found=yes.
+#   my_slu_found
+#      Whether any of the headers.  Either "yes" or "no".
+AC_DEFUN([_slu_search_versions],
+    [
+        my_slu_header=slu_ddefs.h
+        _slu_lib_path($1, $my_slu_header)
+        if test "$my_slu_found" != "yes"; then
+            my_slu_header="dsp_defs.h"
+            _slu_lib_path($1, $my_slu_header)
+        fi
+    ]
+)
+
+
+# _slu_search_default()
+#
+# Search for SuperLU in the default locations "/usr" and "/usr/local".
+#
+# Shell variables:
+#   with_superlu
+#     Root of the SuperLU installation: first of "/usr" and "/usr/local".
+#     Contents is only meaningful for my_slu_found=yes.
+#   For other output variables see documentation of _slu_search_versions().
+AC_DEFUN([_slu_search_default],
+    [
+        with_superlu=/usr
+        _slu_search_versions($with_superlu)
+
+        if test "$my_slu_found" = "no"; then
+            with_superlu=/usr/local
+            _slu_search_versions($with_superlu)
+        fi
+    ]
+)
+
+
+# OPM_PATH_SUPERLU()
+#
+# REQUIRES: AC_PROG_CC, AX_BLAS
+#
+# Shell variables:
+#   with_superlu
+#     "no", "yes (version 4.3 or newer)", "yes (version 4.2 or older, post 2005)" or "yes (pre 2005)"
+#   direct_SUPERLU_CPPFLAGS
+#   direct_SUPERLU_LIBS
+#      CPPFLAGS and LIBS necessary to link against SuperLU.  This variable
+#      contains no indirect references and is suitable for use inside
+#      configure.  Guaranteed empty if SuperLU was not found.
+#   SUPERLU_CPPFLAGS
+#   SUPERLU_LIBS
+#      CPPFLAGS and LIBS necessary to link against SuperLU.  This variable may
+#      contain indirect references and is suitable for use inside makefiles.
+#      Guaranteed empty if SuperLU was not found.
+#   HAVE_SUPERLU
+#      "0" or "1" depending on whether SuperLU was found.
+#
+# Substitutions:
+#   SUPERLU_LIBS
+#   SUPERLU_CPPFLAGS
+#     Substitutes the values of the corresponding shell variables.
+#   ALL_PKG_LIBS
+#   ALL_PKG_CPPFLAGS
+#     Adds references to SuperLU's substitutions.
+#
+# Defines:
+#   HAVE_SUPERLU
+#     ENABLE_SUPERLU or undefined.  Whether SuperLU was found.  The correct
+#     way to check this is "#if HAVE_SUPERLU": This way SuperLU features will
+#     be disabled unless ${SUPERLU_CPPFLAGS} was given when compiling.
+#   SUPERLU_POST_2005_VERSION
+#     1 or undefined.  A post-2005 version of SuperLU uses the header
+#     "slu_ddefs.h" while earlier versions use "dsp_defs.h".
+#   SUPERLU_MIN_VERSION_4_3
+#     1 or undefined. SuperLU version 4.3 or newer uses the symbol
+#     "SLU_DOUBLE" while earlier versions use "DOUBLE".
+#   HAVE_MEM_USAGE_T_EXPANSIONS
+#     1 or undefined.  Whether "mem_usage_t.expansions" was found in
+#     "slu_ddefs.h" or "dsp_defs.h" as apropriate.
+#
+# Conditionals:
+#   SUPERLU
+AC_DEFUN([OPM_PATH_SUPERLU],[
+    AC_REQUIRE([AC_PROG_CC])
+    # we need this for FLIBS
+    AC_REQUIRE([AC_F77_LIBRARY_LDFLAGS])
+    AC_REQUIRE([AX_BLAS])
+
+    #
+    # User hints ...
+    #
+    my_lib_path=""
+    my_include_path=""
+    AC_ARG_WITH([superlu],
+        [AC_HELP_STRING([--with-superlu],[user defined path to SuperLU library])],
+        [dnl
+            if test x"$withval" != xno ; then
+                # get absolute path
+                with_superlu=`eval cd $withval > /dev/null 2>&1 && pwd`
+                if test x"$withval" = xyes; then
+                    # Search in default locations
+                    _slu_search_default
+                else
+                    # Search for the headers in the specified location
+                    _slu_search_versions(["$with_superlu"])
+                fi
+            fi
+        ], [dnl
+            # Search in default locations
+            _slu_search_default
+        ])
+
+    AC_ARG_WITH([superlu-lib],
+        [AC_HELP_STRING([--with-superlu-lib],
+                [The name of the static SuperLU library to link to.  By default
+                    the static library with the name superlu.a is tried, but
+                    only if shared linking has failed first.  Giving this
+                    options forces static linking for SuperLU.])],
+        [
+            if test x"$withval" = xno ; then
+                with_superlu_lib=
+            fi
+        ])
+
+    AC_ARG_WITH([superlu-blaslib],
+        [AC_HELP_STRING([--with-superlu-blaslib],
+                [The name of the static blas library to link to.  By default
+                    we try to link to the systems blas library (see
+                    --with-blas).  Giving this options forces static linking
+                    for SuperLU.])],
+        [
+            if test "$withval" = no ; then
+                with_superlu_blaslib=
+            fi
+        ])
+
+    # store old values
+    ac_save_LDFLAGS="$LDFLAGS"
+    ac_save_CPPFLAGS="$CPPFLAGS"
+    ac_save_LIBS="$LIBS"
+
+    # do nothing if --without-superlu is used
+    if test x"$with_superlu" != x"no" ; then
+        # defaultpath
+        SUPERLU_LIB_PATH="$with_superlu/$my_lib_path"
+        SUPERLU_INCLUDE_PATH="$with_superlu/$my_include_path"
+
+        # set variables so that tests can use them
+        direct_SUPERLU_CPPFLAGS="-I$SUPERLU_INCLUDE_PATH -DENABLE_SUPERLU"
+        SUPERLU_CPPFLAGS="-I$SUPERLU_INCLUDE_PATH -DENABLE_SUPERLU"
+        CPPFLAGS="$CPPFLAGS $direct_SUPERLU_CPPFLAGS"
+
+        # check for central header
+        AC_CHECK_HEADER([$my_slu_header],
+            [HAVE_SUPERLU="1"],
+            [
+                HAVE_SUPERLU="0"
+                AC_MSG_WARN([$my_slu_header not found in $SUPERLU_INCLUDE_PATH with $CPPFLAGS])
+            ])
+
+        # if header is found check for the libs
+        if test x$HAVE_SUPERLU = x1 ; then
+            HAVE_SUPERLU=0
+
+            # if neither --with-superlu-lib nor --with-superlu-blaslib was
+            # given, try to link dynamically or with properly names static libs
+            if test x"$with_superlu_lib$with_superlu_blaslib" = x; then
+                LDFLAGS="$ac_save_LDFLAGS -L$SUPERLU_LIB_PATH"
+                LIBS="$ac_save_LIBS"
+                AC_CHECK_LIB([superlu], [dgssvx],
+                    [
+                        direct_SUPERLU_LIBS="-L$SUPERLU_LIB_PATH -lsuperlu $BLAS_LIBS $FLIBS"
+                        SUPERLU_LIBS="-L$SUPERLU_LIB_PATH -lsuperlu \${BLAS_LIBS} \${FLIBS}"
+                        HAVE_SUPERLU="1"
+                    ], [], [$BLAS_LIBS $FLIBS])
+            fi
+
+            if test $HAVE_SUPERLU = 0 && test x"$with_superlu_lib" = x; then
+                # set the default
+                with_superlu_lib=superlu.a
+            fi
+
+            if test $HAVE_SUPERLU = 0 &&
+                test x"$with_superlu_blaslib" = x; then
+                # try system blas
+                LDFLAGS="$ac_save_LDFLAGS"
+                LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib $BLAS_LIBS $FLIBS $ac_save_LIBS"
+                AC_CHECK_FUNC([dgssvx],
+                    [
+                        direct_SUPERLU_LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib $BLAS_LIBS $FLIBS"
+                        SUPERLU_LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib \${BLAS_LIBS} \${FLIBS}"
+                        HAVE_SUPERLU="1"
+                    ])
+            fi
+
+            # No default for with_superlu_blaslib
+
+            if test $HAVE_SUPERLU = 0 &&
+                    test x"$with_superlu_blaslib" != x; then
+                # try internal blas
+                LDFLAGS="$ac_save_LDFLAGS"
+                LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib $SUPERLU_LIB_PATH/$with_superlu_blaslib $FLIBS $ac_save_LIBS"
+                AC_CHECK_FUNC([dgssvx],
+                    [
+                        direct_SUPERLU_LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib $SUPERLU_LIB_PATH/$with_superlu_blaslib $FLIBS"
+                        SUPERLU_LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib $SUPERLU_LIB_PATH/$with_superlu_blaslib \${FLIBS}"
+                        HAVE_SUPERLU="1"
+                    ])
+            fi
+            
+            # test whether SuperLU version is at least 4.3
+            if test $HAVE_SUPERLU = "1" ; then
+                AC_CHECK_DECL([SLU_DOUBLE], [SUPERLU_MIN_VERSION_4_3="1"], [], [#include <$my_slu_header>])
+            fi
+        fi
+
+    else # $with_superlu = no
+        HAVE_SUPERLU=0
+    fi
+
+    # Inform the user whether SuperLU was sucessfully found
+    AC_MSG_CHECKING([SuperLU])
+    if test x$HAVE_SUPERLU = x1 ; then
+        if test "$my_slu_header" = "slu_ddefs.h"; then
+            if test x$SUPERLU_MIN_VERSION_4_3 = x1 ; then
+                with_superlu="yes (version 4.3 or newer)"
+            else
+                with_superlu="yes (version 4.2 or older, post 2005)"
+            fi
+        else
+            with_superlu="yes (pre 2005)"
+        fi
+    else
+        with_superlu="no"
+    fi
+    AC_MSG_RESULT([$with_superlu])
+
+    # check for optional member
+    if test $HAVE_SUPERLU = 1 ; then
+        if test "$my_slu_header" = "slu_ddefs.h"; then
+            AC_CHECK_MEMBERS([mem_usage_t.expansions],[],[],[#include "slu_ddefs.h"])
+        else
+            AC_CHECK_MEMBERS([mem_usage_t.expansions],[],[],[#include "dsp_defs.h"])
+        fi
+    fi
+
+    # substitute variables
+    if test x$HAVE_SUPERLU = x0 ; then
+        SUPERLU_LIBS=
+        SUPERLU_CPPFLAGS=
+    fi
+    AC_SUBST([SUPERLU_LIBS])
+    AC_SUBST([SUPERLU_CPPFLAGS])
+    #DUNE_ADD_ALL_PKG([SUPERLU], [\${SUPERLU_CPPFLAGS}], [], [\${SUPERLU_LIBS}])
+
+    # tell automake
+    AM_CONDITIONAL(SUPERLU, test x$HAVE_SUPERLU = x1)
+
+    # tell the preprocessor
+    if test x$HAVE_SUPERLU = x1 ; then
+        AC_DEFINE([HAVE_SUPERLU], [ENABLE_SUPERLU], [Define to ENABLE_SUPERLU if SUPERLU is found])
+        if test "$my_slu_header" = "slu_ddefs.h"; then
+            AC_DEFINE([SUPERLU_POST_2005_VERSION], 1, [define to 1 if there is  a header slu_ddefs.h in SuperLU])
+            if test x$SUPERLU_MIN_VERSION_4_3 = x1 ; then
+                AC_DEFINE([SUPERLU_MIN_VERSION_4_3], 1, [define to 1 if there SLU_DOUBLE imported by header slu_ddefs.h from SuperLU])
+            fi
+        fi
+    fi
+
+    # summary
+    #DUNE_ADD_SUMMARY_ENTRY([SuperLU],[$with_superlu])
+
+    # restore variables
+    LDFLAGS="$ac_save_LDFLAGS"
+    CPPFLAGS="$ac_save_CPPFLAGS"
+    LIBS="$ac_save_LIBS"
+])

opm/core/eclipse/EclipseGridParserHelpers.hpp

@@ -398,6 +398,13 @@ namespace
 			yv.push_back(table[k][i]);		    
 		    }
 		}
+                if (xv.empty()) {
+                    // Nothing specified, the entire column is defaulted.
+                    // We insert zeros.
+                    for (int i=0; i<int(indx.size()); ++i) {
+                        table[k][indx[i]] = 0.0;
+                    }
+                } else {
                     // Interpolate
                     for (int i=0; i<int(indx.size()); ++i) {
                         table[k][indx[i]] = linearInterpolationExtrap(xv, yv, x[i]);
@@ -405,6 +412,7 @@ namespace
                 }
 	    }
 	}
+    }
 
     // Reads keywords SGOF and SWOF
     inline void readSGWOF(std::istream& is, table_t& relperm_table,

opm/core/eclipse/SpecialEclipseFields.hpp

@@ -1047,7 +1047,7 @@ struct GCONINJE : public SpecialBase
 	    gconinje_line.injector_type_ = readString(is);
 	    gconinje_line.control_mode_ = readString(is);
 	    std::vector<double> double_data(10, -1.0E20);
-            const int num_to_read = 10;
+            const int num_to_read = 4;
 	    int num_read = readDefaultedVectorData(is, double_data, num_to_read);
 	    gconinje_line.surface_flow_max_rate_ = double_data[0];
             gconinje_line.resv_flow_max_rate_ = double_data[1];

opm/core/fluid/IncompPropertiesInterface.hpp

@@ -62,7 +62,7 @@ namespace Opm
         /// \return Array of P viscosity values.
         virtual const double* viscosity() const = 0;
 
-        /// Densities of fluid phases at surface conditions.
+        /// Densities of fluid phases at reservoir conditions.
         /// \return Array of P density values.
         virtual const double* density() const = 0;
 

opm/core/grid/cpgpreprocess/preprocess.c

@@ -590,6 +590,9 @@ get_zcorn_sign(int nx, int ny, int nz, const int *actnum,
                     c1 = i/2 + nx*(j/2 + ny*(k/2));
                     c2 = i/2 + nx*(j/2 + ny*((k+1)/2));
 
+                    assert (c1 < (nx * ny * nz));
+                    assert (c2 < (nx * ny * nz));
+
                     if (((actnum == NULL) ||
                          (actnum[c1] && actnum[c2]))
                         && (z2 < z1)) {

opm/core/newwells.c

@@ -468,6 +468,8 @@ add_well(enum WellType  type     ,
 
     struct WellMgmt *m;
 
+    assert (W != NULL);
+
     nw        = W->number_of_wells;
     nperf_tot = W->well_connpos[nw];
 
@@ -532,6 +534,7 @@ append_well_controls(enum WellControlType type,
     struct WellControlMgmt *m;
 
     assert (W != NULL);
+    assert ((0 <= well_index) && (well_index < W->number_of_wells));
 
     ctrl = W->ctrls[well_index];
     np = W->number_of_phases;
@@ -567,8 +570,13 @@ void
 set_current_control(int well_index, int current_control, struct Wells *W)
 /* ---------------------------------------------------------------------- */
 {
+    assert (W != NULL);
+    assert ((0 <= well_index) && (well_index < W->number_of_wells));
+
     assert (W->ctrls[well_index] != NULL);
 
+    assert (current_control < W->ctrls[well_index]->num);
+
     W->ctrls[well_index]->current = current_control;
 }
 
@@ -578,7 +586,85 @@ void
 clear_well_controls(int well_index, struct Wells *W)
 /* ---------------------------------------------------------------------- */
 {
+    assert (W != NULL);
+    assert ((0 <= well_index) && (well_index < W->number_of_wells));
+
     if (W->ctrls[well_index] != NULL) {
         W->ctrls[well_index]->num = 0;
     }
 }
+
+
+/* ---------------------------------------------------------------------- */
+struct Wells *
+clone_wells(const struct Wells *W)
+/* ---------------------------------------------------------------------- */
+{
+    int                        c, np, nperf, ok, pos, w;
+    double                     target;
+    const int                 *cells;
+    const double              *WI, *comp_frac, *distr;
+    enum WellControlType       type;
+    const struct WellControls *ctrls;
+
+    struct Wells *ret;
+
+    if (W == NULL) {
+        ret = NULL;
+    }
+    else {
+        np  = W->number_of_phases;
+        ret = create_wells(W->number_of_phases, W->number_of_wells,
+                           W->well_connpos[ W->number_of_wells ]);
+
+        if (ret != NULL) {
+            pos = W->well_connpos[ 0 ];
+            ok  = 1;
+
+            for (w = 0; ok && (w < W->number_of_wells); w++) {
+                nperf = W->well_connpos[w + 1] - pos;
+                cells = W->well_cells + pos;
+
+                WI        = W->WI        != NULL ? W->WI        + pos  : NULL;
+                comp_frac = W->comp_frac != NULL ? W->comp_frac + w*np : NULL;
+
+                ok = add_well(W->type[ w ], W->depth_ref[ w ], nperf,
+                              comp_frac, cells, WI, W->name[ w ], ret);
+
+                /* Capacity should be sufficient from create_wells() */
+                assert (ok);
+
+                ctrls = W->ctrls[ w ];
+
+                if (ok) {
+                    ok = well_controls_reserve(ctrls->num, np, ret->ctrls[ w ]);
+
+                    for (c = 0; ok && (c < ctrls->num); c++) {
+                        type   =   ctrls->type  [ c ];
+                        target =   ctrls->target[ c ];
+                        distr  = & ctrls->distr [ c * np ];
+
+                        ok = append_well_controls(type, target, distr, w, ret);
+
+                        /* Capacity *should* be sufficient from
+                         *  well_controls_reserve() */
+                        assert (ok);
+                    }
+                }
+
+                if (ok) {
+                    set_current_control(w, ctrls->current, ret);
+                }
+
+                pos = W->well_connpos[w + 1];
+            }
+
+            if (! ok) {
+                destroy_wells(ret);
+                ret = NULL;
+            }
+        }
+    }
+
+    return ret;
+}

opm/core/newwells.h

@@ -235,6 +235,19 @@ void
 destroy_wells(struct Wells *W);
 
 
+/**
+ * Create a deep-copy (i.e., clone) of an existing Wells object, including its
+ * controls.
+ *
+ * @param[in] W Existing Wells object.
+ * @return Complete clone of the input object.  Dispose of resources using
+ * function destroy_wells() when no longer needed.  Returns @c NULL in case of
+ * allocation failure.
+ */
+struct Wells *
+clone_wells(const struct Wells *W);
+
+
 #ifdef __cplusplus
 }
 #endif

opm/core/simulator/SimulatorCompressibleTwophase.cpp

@@ -407,8 +407,8 @@ namespace Opm
                 // Optionally, check if well controls are satisfied.
                 if (check_well_controls_) {
                     Opm::computePhaseFlowRatesPerWell(*wells_,
-                                                      fractional_flows,
                                                       well_state.perfRates(),
+                                                      fractional_flows,
                                                       well_resflows_phase);
                     std::cout << "Checking well conditions." << std::endl;
                     // For testing we set surface := reservoir

opm/core/simulator/SimulatorIncompTwophase.cpp

@@ -485,8 +485,8 @@ namespace Opm
                 // Optionally, check if well controls are satisfied.
                 if (check_well_controls_) {
                     Opm::computePhaseFlowRatesPerWell(*wells_,
-                                                      fractional_flows,
                                                       well_state.perfRates(),
+                                                      fractional_flows,
                                                       well_resflows_phase);
                     std::cout << "Checking well conditions." << std::endl;
                     // For testing we set surface := reservoir

opm/core/transport/reorder/TransportModelInterface.cpp

@@ -20,6 +20,7 @@
 #include <opm/core/transport/reorder/TransportModelInterface.hpp>
 #include <opm/core/transport/reorder/reordersequence.h>
 #include <opm/core/grid.h>
+#include <opm/core/utility/StopWatch.hpp>
 
 #include <vector>
 #include <cassert>
@@ -31,7 +32,11 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g
     std::vector<int> sequence(grid.number_of_cells);
     std::vector<int> components(grid.number_of_cells + 1);
     int ncomponents;
+    time::StopWatch clock;
+    clock.start();
     compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents);
+    clock.stop();
+    std::cout << "Topological sort took: " << clock.secsSinceStart() << " seconds." << std::endl;
 
     // Invoke appropriate solve method for each interdependent component.
     for (int comp = 0; comp < ncomponents; ++comp) {

opm/core/transport/reorder/TransportModelTracerTof.cpp

@@ -0,0 +1,122 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/core/transport/reorder/TransportModelTracerTof.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <algorithm>
+#include <numeric>
+#include <cmath>
+
+namespace Opm
+{
+
+
+    /// Construct solver.
+    /// \param[in] grid      A 2d or 3d grid.
+    TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid)
+        : grid_(grid)
+    {
+    }
+
+
+
+
+    /// Solve for time-of-flight.
+    /// \param[in]  darcyflux         Array of signed face fluxes.
+    /// \param[in]  porevolume        Array of pore volumes.
+    /// \param[in]  source            Source term. Sign convention is:
+    ///                                 (+) inflow flux,
+    ///                                 (-) outflow flux.
+    /// \param[out] tof               Array of time-of-flight values.
+    void TransportModelTracerTof::solveTof(const double* darcyflux,
+                                           const double* porevolume,
+                                           const double* source,
+                                           std::vector<double>& tof)
+    {
+        darcyflux_ = darcyflux;
+        porevolume_ = porevolume;
+        source_ = source;
+#ifndef NDEBUG
+        // Sanity check for sources.
+        const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
+        if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
+            THROW("Sources do not sum to zero: " << cum_src);
+        }
+#endif
+        tof.resize(grid_.number_of_cells);
+        std::fill(tof.begin(), tof.end(), 0.0);
+        tof_ = &tof[0];
+        reorderAndTransport(grid_, darcyflux);
+    }
+
+
+
+
+    void TransportModelTracerTof::solveSingleCell(const int cell)
+    {
+        // Compute flux terms.
+        // Sources have zero tof, and therefore do not contribute
+        // to upwind_term. Sinks on the other hand, must be added
+        // to the downwind_flux (note sign change resulting from
+        // different sign conventions: pos. source is injection,
+        // pos. flux is outflow).
+        double upwind_term = 0.0;
+        double downwind_flux = std::max(-source_[cell], 0.0);
+        for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
+            int f = grid_.cell_faces[i];
+            double flux;
+            int other;
+            // Compute cell flux
+            if (cell == grid_.face_cells[2*f]) {
+                flux  = darcyflux_[f];
+                other = grid_.face_cells[2*f+1];
+            } else {
+                flux  =-darcyflux_[f];
+                other = grid_.face_cells[2*f];
+            }
+            // Add flux to upwind_term or downwind_flux, if interior.
+            if (other != -1) {
+                if (flux < 0.0) {
+                    upwind_term  += flux*tof_[other];
+                } else {
+                    downwind_flux += flux;
+                }
+            }
+        }
+
+        // Compute tof.
+        tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux;
+    }
+
+
+
+
+    void TransportModelTracerTof::solveMultiCell(const int num_cells, const int* cells)
+    {
+        std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
+        for (int i = 0; i < num_cells; ++i) {
+            solveSingleCell(cells[i]);
+        }
+    }
+
+
+
+
+} // namespace Opm

opm/core/transport/reorder/TransportModelTracerTof.hpp

@@ -0,0 +1,74 @@
+/*
+  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_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED
+#define OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED
+
+#include <opm/core/transport/reorder/TransportModelInterface.hpp>
+#include <vector>
+#include <map>
+#include <ostream>
+struct UnstructuredGrid;
+
+namespace Opm
+{
+
+    class IncompPropertiesInterface;
+
+    /// Implements a first-order finite volume solver for
+    /// (single-phase) time-of-flight using reordering.
+    /// The equation solved is:
+    ///     v \cdot \grad\tau = \phi
+    /// where v is the fluid velocity, \tau is time-of-flight and
+    /// \phi is the porosity. This is a boundary value problem, where
+    /// \tau is specified to be zero on all inflow boundaries.
+    class TransportModelTracerTof : public TransportModelInterface
+    {
+    public:
+        /// Construct solver.
+        /// \param[in] grid      A 2d or 3d grid.
+        TransportModelTracerTof(const UnstructuredGrid& grid);
+
+        /// Solve for time-of-flight.
+        /// \param[in]  darcyflux         Array of signed face fluxes.
+        /// \param[in]  porevolume        Array of pore volumes.
+        /// \param[in]  source            Source term. Sign convention is:
+        ///                                 (+) inflow flux,
+        ///                                 (-) outflow flux.
+        /// \param[out] tof               Array of time-of-flight values.
+        void solveTof(const double* darcyflux,
+                      const double* porevolume,
+                      const double* source,
+                      std::vector<double>& tof);
+
+    private:
+        virtual void solveSingleCell(const int cell);
+        virtual void solveMultiCell(const int num_cells, const int* cells);
+
+    private:
+        const UnstructuredGrid& grid_;
+        const double* darcyflux_;   // one flux per grid face
+        const double* porevolume_;  // one volume per cell
+        const double* source_;      // one volumetric source term per cell
+        double* tof_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED

opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp

@@ -0,0 +1,671 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/linalg/blas_lapack.h>
+#include <algorithm>
+#include <numeric>
+#include <cmath>
+
+namespace Opm
+{
+
+
+    // ---------------   Helpers for TransportModelTracerTofDiscGal ---------------
+
+
+
+    /// A class providing discontinuous Galerkin basis functions.
+    struct DGBasis
+    {
+        static int numBasisFunc(const int dimensions,
+                                const int degree)
+        {
+            switch (dimensions) {
+            case 1:
+                return degree + 1;
+            case 2:
+                return (degree + 2)*(degree + 1)/2;
+            case 3:
+                return (degree + 3)*(degree + 2)*(degree + 1)/6;
+            default:
+                THROW("Dimensions must be 1, 2 or 3.");
+            }
+        }
+
+        /// Evaluate all nonzero basis functions at x,
+        /// writing to f_x. The array f_x must have
+        /// size numBasisFunc(grid.dimensions, degree).
+        ///
+        /// The basis functions are the following
+        ///     Degree 0: 1.
+        ///     Degree 1: x - xc, y - yc, z - zc etc.
+        /// Further degrees await development.
+        static void eval(const UnstructuredGrid& grid,
+                         const int cell,
+                         const int degree,
+                         const double* x,
+                         double* f_x)
+        {
+            const int dim = grid.dimensions;
+            const double* cc = grid.cell_centroids + dim*cell;
+            // Note intentional fallthrough in this switch statement!
+            switch (degree) {
+            case 1:
+                for (int ix = 0; ix < dim; ++ix) {
+                    f_x[1 + ix] = x[ix] - cc[ix];
+                }
+            case 0:
+                f_x[0] = 1;
+                break;
+            default:
+                THROW("Maximum degree is 1 for now.");
+            }
+        }
+
+        /// Evaluate gradients of all nonzero basis functions at x,
+        /// writing to grad_f_x. The array grad_f_x must have size
+        /// numBasisFunc(grid.dimensions, degree) * grid.dimensions.
+        /// The <grid.dimensions> components of the first basis function
+        /// gradient come before the components of the second etc.
+        static void evalGrad(const UnstructuredGrid& grid,
+                             const int /*cell*/,
+                             const int degree,
+                             const double* /*x*/,
+                             double* grad_f_x)
+        {
+            const int dim = grid.dimensions;
+            const int num_basis = numBasisFunc(dim, degree);
+            std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0);
+            if (degree > 1) {
+                THROW("Maximum degree is 1 for now.");
+            } else if (degree == 1) {
+                for (int ix = 0; ix < dim; ++ix) {
+                    grad_f_x[dim*(ix + 1) + ix] = 1.0;
+                }
+            }
+        }
+    };
+
+
+
+    static void cross(const double* a, const double* b, double* res)
+    {
+        res[0] = a[1]*b[2] - a[2]*b[1];
+        res[1] = a[2]*b[0] - a[0]*b[2];
+        res[2] = a[0]*b[1] - a[1]*b[0];
+    }
+
+
+
+
+    static double triangleArea3d(const double* p0,
+                                 const double* p1,
+                                 const double* p2)
+    {
+        double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
+        double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
+        double cr[3];
+        cross(a, b, cr);
+        return 0.5*std::sqrt(cr[0]*cr[0] + cr[1]*cr[1] + cr[2]*cr[2]);
+    }
+
+
+
+
+    /// Calculates the determinant of a 3 x 3 matrix, represented as
+    /// three three-dimensional arrays.
+    static double determinantOf(const double* a0,
+                                const double* a1,
+                                const double* a2)
+    {
+        return
+            a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) -
+            a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) +
+            a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]);
+    }
+
+
+
+
+    /// Computes the volume of a tetrahedron consisting of 4 vertices
+    /// with 3-dimensional coordinates
+    static double tetVolume(const double* p0,
+                            const double* p1,
+                            const double* p2,
+                            const double* p3)
+    {
+        double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
+        double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
+        double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] };
+        return std::fabs(determinantOf(a, b, c) / 6.0);
+    }
+
+
+
+
+    /// A class providing numerical quadrature for cells.
+    /// In general: \int_{cell} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
+    /// Note that this class does multiply weights by cell volume,
+    /// so weights always sum to cell volume.
+    /// Degree 1 method:
+    ///     Midpoint (centroid) method.
+    ///         n = 1, w_0 = cell volume, x_0 = cell centroid
+    /// Degree 2 method:
+    ///    Based on subdivision of each cell face into triangles
+    ///    with the face centroid as a common vertex, and then
+    ///    subdividing the cell into tetrahedra with the cell
+    ///    centroid as a common vertex. Then apply the tetrahedron
+    ///    rule with the following 4 nodes (uniform weights):
+    ///        a = 0.138196601125010515179541316563436
+    ///        x_i has all barycentric coordinates = a, except for
+    ///            the i'th coordinate which is = 1 - 3a.
+    ///    This rule is from http://nines.cs.kuleuven.be/ecf,
+    ///    it is the second degree 2 4-point rule for tets,
+    ///    referenced to Stroud(1971).
+    ///    The tetrahedra are numbered T_{i,j}, and are given by the
+    ///    cell centroid C, the face centroid FC_i, and two nodes
+    ///    of face i: FN_{i,j}, FN_{i,j+1}.
+    class CellQuadrature
+    {
+    public:
+        CellQuadrature(const UnstructuredGrid& grid,
+                       const int cell,
+                       const int degree)
+            : grid_(grid), cell_(cell), degree_(degree)
+        {
+            if (degree > 2) {
+                THROW("CellQuadrature exact for polynomial degrees > 1 not implemented.");
+            }
+            if (degree == 2) {
+                // Prepare subdivision.
+            }
+        }
+
+        int numQuadPts() const
+        {
+            if (degree_ < 2) {
+                return 1;
+            }
+            // Degree 2 case.
+            int sumnodes = 0;
+            for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
+                const int face = grid_.cell_faces[hf];
+                sumnodes += grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
+            }
+            return 4*sumnodes;
+        }
+
+        void quadPtCoord(const int index, double* coord) const
+        {
+            const int dim = grid_.dimensions;
+            const double* cc = grid_.cell_centroids + dim*cell_;
+            if (degree_ < 2) {
+                std::copy(cc, cc + dim, coord);
+                return;
+            }
+            // Degree 2 case.
+            int tetindex = index / 4;
+            const int subindex = index % 4;
+            const double* nc = grid_.node_coordinates;
+            for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
+                const int face = grid_.cell_faces[hf];
+                const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
+                if (nfn <= tetindex) {
+                    // Our tet is not associated with this face.
+                    tetindex -= nfn;
+                    continue;
+                }
+                const double* fc = grid_.face_centroids + dim*face;
+                const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
+                const int node0 = fnodes[tetindex];
+                const int node1 = fnodes[(tetindex + 1) % nfn];
+                const double* n0c = nc + dim*node0;
+                const double* n1c = nc + dim*node1;
+                const double a = 0.138196601125010515179541316563436;
+                // Barycentric coordinates of our point in the tet.
+                double baryc[4] = { a, a, a, a };
+                baryc[subindex] = 1.0 - 3.0*a;
+                for (int dd = 0; dd < dim; ++dd) {
+                    coord[dd] = baryc[0]*cc[dd] + baryc[1]*fc[dd] + baryc[2]*n0c[dd] + baryc[3]*n1c[dd];
+                }
+                return;
+            }
+            THROW("Should never reach this point.");
+        }
+
+        double quadPtWeight(const int index) const
+        {
+            if (degree_ < 2) {
+                return grid_.cell_volumes[cell_];
+            }
+            // Degree 2 case.
+            const int dim = grid_.dimensions;
+            const double* cc = grid_.cell_centroids + dim*cell_;
+            int tetindex = index / 4;
+            const double* nc = grid_.node_coordinates;
+            for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
+                const int face = grid_.cell_faces[hf];
+                const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
+                if (nfn <= tetindex) {
+                    // Our tet is not associated with this face.
+                    tetindex -= nfn;
+                    continue;
+                }
+                const double* fc = grid_.face_centroids + dim*face;
+                const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
+                const int node0 = fnodes[tetindex];
+                const int node1 = fnodes[(tetindex + 1) % nfn];
+                const double* n0c = nc + dim*node0;
+                const double* n1c = nc + dim*node1;
+                return 0.25*tetVolume(cc, fc, n0c, n1c);
+            }
+            THROW("Should never reach this point.");
+        }
+
+    private:
+        const UnstructuredGrid& grid_;
+        const int cell_;
+        const int degree_;
+    };
+
+
+
+
+
+    /// A class providing numerical quadrature for faces.
+    /// In general: \int_{face} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
+    /// Note that this class does multiply weights by face area,
+    /// so weights always sum to face area.
+    /// Degree 1 method:
+    ///     Midpoint (centroid) method.
+    ///         n = 1, w_0 = face area, x_0 = face centroid
+    /// Degree 2 method:
+    ///    Based on subdivision of the face into triangles,
+    ///    with the centroid as a common vertex, and the triangle
+    ///    edge midpoint rule.
+    ///    Triangle i consists of the centroid C, nodes N_i and N_{i+1}.
+    ///    Its area is A_i.
+    ///        n = 2 * nn  (nn = num nodes in face)
+    ///        For i = 0..(nn-1):
+    ///        w_i      = 1/3 A_i.
+    ///        w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i
+    ///        x_i      = (N_i + N_{i+1})/2
+    ///        x_{nn+i} = (C + N_i)/2
+    ///    All N and A indices are interpreted cyclic, modulus nn.
+    class FaceQuadrature
+    {
+    public:
+        FaceQuadrature(const UnstructuredGrid& grid,
+                       const int face,
+                       const int degree)
+            : grid_(grid), face_(face), degree_(degree)
+        {
+            if (grid_.dimensions != 3) {
+                THROW("FaceQuadrature only implemented for 3D case.");
+            }
+            if (degree_ > 2) {
+                THROW("FaceQuadrature exact for polynomial degrees > 2 not implemented.");
+            }
+        }
+
+        int numQuadPts() const
+        {
+            if (degree_ < 2) {
+                return 1;
+            }
+            // Degree 2 case.
+            return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]);
+        }
+
+        void quadPtCoord(const int index, double* coord) const
+        {
+            const int dim = grid_.dimensions;
+            const double* fc = grid_.face_centroids + dim*face_;
+            if (degree_ < 2) {
+                std::copy(fc, fc + dim, coord);
+                return;
+            }
+            // Degree 2 case.
+            const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
+            const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
+            const double* nc = grid_.node_coordinates;
+            if (index < nn) {
+                // Boundary edge midpoint.
+                const int node0 = fnodes[index];
+                const int node1 = fnodes[(index + 1)%nn];
+                for (int dd = 0; dd < dim; ++dd) {
+                    coord[dd] = 0.5*(nc[dim*node0 + dd] + nc[dim*node1 + dd]);
+                }
+            } else {
+                // Interiour edge midpoint.
+                // Recall that index is now in [nn, 2*nn).
+                const int node = fnodes[index - nn];
+                for (int dd = 0; dd < dim; ++dd) {
+                    coord[dd] = 0.5*(nc[dim*node + dd] + fc[dd]);
+                }
+            }
+        }
+
+        double quadPtWeight(const int index) const
+        {
+            if (degree_ < 2) {
+                return grid_.face_areas[face_];
+            }
+            // Degree 2 case.
+            const int dim = grid_.dimensions;
+            const double* fc = grid_.face_centroids + dim*face_;
+            const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
+            const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
+            const double* nc = grid_.node_coordinates;
+            if (index < nn) {
+                // Boundary edge midpoint.
+                const int node0 = fnodes[index];
+                const int node1 = fnodes[(index + 1)%nn];
+                const double area = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
+                return area / 3.0;
+            } else {
+                // Interiour edge midpoint.
+                // Recall that index is now in [nn, 2*nn).
+                const int node0 = fnodes[(index - 1) % nn];
+                const int node1 = fnodes[index - nn];
+                const int node2 = fnodes[(index + 1) % nn];
+                const double area0 = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
+                const double area1 = triangleArea3d(nc + dim*node2, nc + dim*node1, fc);
+                return (area0 + area1) / 3.0;
+            }
+        }
+
+    private:
+        const UnstructuredGrid& grid_;
+        const int face_;
+        const int degree_;
+    };
+
+
+
+
+    // Initial version: only a constant interpolation.
+    static void interpolateVelocity(const UnstructuredGrid& grid,
+                                    const int cell,
+                                    const double* darcyflux,
+                                    const double* /*x*/,
+                                    double* v)
+    {
+        const int dim = grid.dimensions;
+        std::fill(v, v + dim, 0.0);
+        const double* cc = grid.cell_centroids + cell*dim;
+        for (int hface = grid.cell_facepos[cell]; hface < grid.cell_facepos[cell+1]; ++hface) {
+            const int face = grid.cell_faces[hface];
+            const double* fc = grid.face_centroids + face*dim;
+            double flux = 0.0;
+            if (cell == grid.face_cells[2*face]) {
+                flux = darcyflux[face];
+            } else {
+                flux = -darcyflux[face];
+            }
+            for (int dd = 0; dd < dim; ++dd) {
+                v[dd] += flux * (fc[dd] - cc[dd]) / grid.cell_volumes[cell];
+            }
+        }
+    }
+
+
+
+    // ---------------   Methods of TransportModelTracerTofDiscGal ---------------
+
+
+
+    /// Construct solver.
+    /// \param[in] grid      A 2d or 3d grid.
+    TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid)
+        : grid_(grid),
+          coord_(grid.dimensions),
+          velocity_(grid.dimensions)
+    {
+    }
+
+
+
+
+    /// Solve for time-of-flight.
+    /// \param[in]  darcyflux         Array of signed face fluxes.
+    /// \param[in]  porevolume        Array of pore volumes.
+    /// \param[in]  source            Source term. Sign convention is:
+    ///                                 (+) inflow flux,
+    ///                                 (-) outflow flux.
+    /// \param[in]  degree            Polynomial degree of DG basis functions used.
+    /// \param[out] tof_coeff         Array of time-of-flight solution coefficients.
+    ///                               The values are ordered by cell, meaning that
+    ///                               the K coefficients corresponding to the first
+    ///                               cell comes before the K coefficients corresponding
+    ///                               to the second cell etc.
+    ///                               K depends on degree and grid dimension.
+    void TransportModelTracerTofDiscGal::solveTof(const double* darcyflux,
+                                                  const double* porevolume,
+                                                  const double* source,
+                                                  const int degree,
+                                                  std::vector<double>& tof_coeff)
+    {
+        darcyflux_ = darcyflux;
+        porevolume_ = porevolume;
+        source_ = source;
+#ifndef NDEBUG
+        // Sanity check for sources.
+        const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
+        if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
+            THROW("Sources do not sum to zero: " << cum_src);
+        }
+#endif
+        degree_ = degree;
+        const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_);
+        tof_coeff.resize(num_basis*grid_.number_of_cells);
+        std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0);
+        tof_coeff_ = &tof_coeff[0];
+        rhs_.resize(num_basis);
+        jac_.resize(num_basis*num_basis);
+        basis_.resize(num_basis);
+        basis_nb_.resize(num_basis);
+        grad_basis_.resize(num_basis*grid_.dimensions);
+        reorderAndTransport(grid_, darcyflux);
+    }
+
+
+
+
+    void TransportModelTracerTofDiscGal::solveSingleCell(const int cell)
+    {
+        // Residual:
+        // For each cell K, basis function b_j (spanning V_h),
+        // writing the solution u_h|K = \sum_i c_i b_i
+        //  Res = - \int_K \sum_i c_i b_i v(x) \cdot \grad b_j dx
+        //        + \int_{\partial K} F(u_h, u_h^{ext}, v(x) \cdot n) b_j ds
+        //        - \int_K \phi b_j
+        // This is linear in c_i, so we do not need any nonlinear iterations.
+        // We assemble the jacobian and the right-hand side. The residual is
+        // equal to Res = Jac*c - rhs, and we compute rhs directly.
+
+        const int dim = grid_.dimensions;
+        const int num_basis = DGBasis::numBasisFunc(dim, degree_);
+
+        std::fill(rhs_.begin(), rhs_.end(), 0.0);
+        std::fill(jac_.begin(), jac_.end(), 0.0);
+
+        // Compute cell residual contribution.
+        // Note: Assumes that \int_K b_j = 0 for all j > 0
+        rhs_[0] += porevolume_[cell];
+
+        // Compute upstream residual contribution.
+        for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
+            const int face = grid_.cell_faces[hface];
+            double flux = 0.0;
+            int upstream_cell = -1;
+            if (cell == grid_.face_cells[2*face]) {
+                flux = darcyflux_[face];
+                upstream_cell = grid_.face_cells[2*face+1];
+            } else {
+                flux = -darcyflux_[face];
+                upstream_cell = grid_.face_cells[2*face];
+            }
+            if (upstream_cell < 0) {
+                // This is an outer boundary. Assumed tof = 0 on inflow, so no contribution.
+                continue;
+            }
+            if (flux >= 0.0) {
+                // This is an outflow boundary.
+                continue;
+            }
+            // Do quadrature over the face to compute
+            // \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds
+            // (where u_h^{ext} is the upstream unknown (tof)).
+            const double normal_velocity = flux / grid_.face_areas[face];
+            FaceQuadrature quad(grid_, face, degree_);
+            for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
+                quad.quadPtCoord(quad_pt, &coord_[0]);
+                DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
+                DGBasis::eval(grid_, upstream_cell, degree_, &coord_[0], &basis_nb_[0]);
+                const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(),
+                                                               tof_coeff_ + num_basis*upstream_cell, 0.0);
+                const double w = quad.quadPtWeight(quad_pt);
+                for (int j = 0; j < num_basis; ++j) {
+                    rhs_[j] -= w * tof_upstream * normal_velocity * basis_[j];
+                }
+            }
+        }
+
+        // Compute cell jacobian contribution. We use Fortran ordering
+        // for jac_, i.e. rows cycling fastest.
+        {
+            CellQuadrature quad(grid_, cell, 2*degree_ - 1);
+            for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
+                // b_i (v \cdot \grad b_j)
+                quad.quadPtCoord(quad_pt, &coord_[0]);
+                DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
+                DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]);
+                interpolateVelocity(grid_, cell, darcyflux_, &coord_[0], &velocity_[0]);
+                const double w = quad.quadPtWeight(quad_pt);
+                for (int j = 0; j < num_basis; ++j) {
+                    for (int i = 0; i < num_basis; ++i) {
+                        for (int dd = 0; dd < dim; ++dd) {
+                            jac_[j*num_basis + i] -= w * basis_[j] * grad_basis_[dim*i + dd] * velocity_[dd];
+                        }
+                    }
+                }
+            }
+        }
+
+        // Compute downstream jacobian contribution from faces.
+        for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
+            const int face = grid_.cell_faces[hface];
+            double flux = 0.0;
+            if (cell == grid_.face_cells[2*face]) {
+                flux = darcyflux_[face];
+            } else {
+                flux = -darcyflux_[face];
+            }
+            if (flux <= 0.0) {
+                // This is an inflow boundary.
+                continue;
+            }
+            // Do quadrature over the face to compute
+            // \int_{\partial K} b_i (v(x) \cdot n) b_j ds
+            const double normal_velocity = flux / grid_.face_areas[face];
+            FaceQuadrature quad(grid_, face, 2*degree_);
+            for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
+                // u^ext flux B   (B = {b_j})
+                quad.quadPtCoord(quad_pt, &coord_[0]);
+                DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
+                const double w = quad.quadPtWeight(quad_pt);
+                for (int j = 0; j < num_basis; ++j) {
+                    for (int i = 0; i < num_basis; ++i) {
+                        jac_[j*num_basis + i] += w * basis_[i] * normal_velocity * basis_[j];
+                    }
+                }
+            }
+        }
+
+        // Compute downstream jacobian contribution from sink terms.
+        // Contribution from inflow sources would be
+        // similar to the contribution from upstream faces, but
+        // it is zero since we let all external inflow be associated
+        // with a zero tof.
+        if (source_[cell] < 0.0) {
+            // A sink.
+            const double flux = -source_[cell]; // Sign convention for flux: outflux > 0.
+            const double flux_density = flux / grid_.cell_volumes[cell];
+            // Do quadrature over the cell to compute
+            // \int_{K} b_i flux b_j dx
+            CellQuadrature quad(grid_, cell, 2*degree_);
+            for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
+                quad.quadPtCoord(quad_pt, &coord_[0]);
+                DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
+                const double w = quad.quadPtWeight(quad_pt);
+                for (int j = 0; j < num_basis; ++j) {
+                    for (int i = 0; i < num_basis; ++i) {
+                        jac_[j*num_basis + i] += w * basis_[i] * flux_density * basis_[j];
+                    }
+                }
+            }
+        }
+
+        // Solve linear equation.
+        MAT_SIZE_T n = num_basis;
+        MAT_SIZE_T nrhs = 1;
+        MAT_SIZE_T lda = num_basis;
+        std::vector<MAT_SIZE_T> piv(num_basis);
+        MAT_SIZE_T ldb = num_basis;
+        MAT_SIZE_T info = 0;
+        dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info);
+        if (info != 0) {
+            // Print the local matrix and rhs.
+            std::cerr << "Failed solving single-cell system Ax = b in cell " << cell
+                      << " with A = \n";
+            for (int row = 0; row < n; ++row) {
+                for (int col = 0; col < n; ++col) {
+                    std::cerr << "    " << jac_[row + n*col];
+                }
+                std::cerr << '\n';
+            }
+            std::cerr << "and b = \n";
+            for (int row = 0; row < n; ++row) {
+                std::cerr << "    " << rhs_[row] << '\n';
+            }
+            THROW("Lapack error: " << info << " encountered in cell " << cell);
+        }
+        // The solution ends up in rhs_, so we must copy it.
+        std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell);
+    }
+
+
+
+
+    void TransportModelTracerTofDiscGal::solveMultiCell(const int num_cells, const int* cells)
+    {
+        std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
+        for (int i = 0; i < num_cells; ++i) {
+            solveSingleCell(cells[i]);
+        }
+    }
+
+
+
+
+} // namespace Opm

opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp

@@ -0,0 +1,93 @@
+/*
+  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_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
+#define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
+
+#include <opm/core/transport/reorder/TransportModelInterface.hpp>
+#include <vector>
+#include <map>
+#include <ostream>
+struct UnstructuredGrid;
+
+namespace Opm
+{
+
+    class IncompPropertiesInterface;
+
+    /// Implements a discontinuous Galerkin solver for
+    /// (single-phase) time-of-flight using reordering.
+    /// The equation solved is:
+    ///     v \cdot \grad\tau = \phi
+    /// where v is the fluid velocity, \tau is time-of-flight and
+    /// \phi is the porosity. This is a boundary value problem, where
+    /// \tau is specified to be zero on all inflow boundaries.
+    /// The user may specify the polynomial degree of the basis function space
+    /// used, but only degrees 0 and 1 are supported so far.
+    class TransportModelTracerTofDiscGal : public TransportModelInterface
+    {
+    public:
+        /// Construct solver.
+        /// \param[in] grid      A 2d or 3d grid.
+        TransportModelTracerTofDiscGal(const UnstructuredGrid& grid);
+
+
+        /// Solve for time-of-flight.
+        /// \param[in]  darcyflux         Array of signed face fluxes.
+        /// \param[in]  porevolume        Array of pore volumes.
+        /// \param[in]  source            Source term. Sign convention is:
+        ///                                 (+) inflow flux,
+        ///                                 (-) outflow flux.
+        /// \param[in]  degree            Polynomial degree of DG basis functions used.
+        /// \param[out] tof_coeff         Array of time-of-flight solution coefficients.
+        ///                               The values are ordered by cell, meaning that
+        ///                               the K coefficients corresponding to the first
+        ///                               cell comes before the K coefficients corresponding
+        ///                               to the second cell etc.
+        ///                               K depends on degree and grid dimension.
+        void solveTof(const double* darcyflux,
+                      const double* porevolume,
+                      const double* source,
+                      const int degree,
+                      std::vector<double>& tof_coeff);
+
+    private:
+        virtual void solveSingleCell(const int cell);
+        virtual void solveMultiCell(const int num_cells, const int* cells);
+
+    private:
+        const UnstructuredGrid& grid_;
+        const double* darcyflux_;   // one flux per grid face
+        const double* porevolume_;  // one volume per cell
+        const double* source_;      // one volumetric source term per cell
+        int degree_;
+        double* tof_coeff_;
+        std::vector<double> rhs_;   // single-cell right-hand-side
+        std::vector<double> jac_;   // single-cell jacobian
+        // Below: storage for quantities needed by solveSingleCell().
+        std::vector<double> coord_;
+        std::vector<double> basis_;
+        std::vector<double> basis_nb_;
+        std::vector<double> grad_basis_;
+        std::vector<double> velocity_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED

opm/core/utility/miscUtilities.cpp

@@ -553,6 +553,7 @@ namespace Opm
     {
         const int np = wells.number_of_phases;
         const int nw = wells.number_of_wells;
+        ASSERT(int(flow_rates_per_well_cell.size()) == wells.well_connpos[nw]);
         phase_flow_per_well.resize(nw * np);
         for (int wix = 0; wix < nw; ++wix) {
             for (int phase = 0; phase < np; ++phase) {

opm/core/wells/WellsGroup.cpp

@@ -572,6 +572,7 @@ namespace Opm
                 break;
             }
             const double total_produced = getTotalProductionFlow(well_surfacerates_phase, phase);
+            const double total_reinjected = - total_produced; // Production negative, injection positive
             const double my_guide_rate = injectionGuideRate(true);
             for (size_t i = 0; i < children_.size(); ++i) {
                 // Apply for all children. 
@@ -580,11 +581,11 @@ namespace Opm
                 const double children_guide_rate = children_[i]->injectionGuideRate(true);
 #ifdef DIRTY_WELLCTRL_HACK
                 children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
-                        (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
+                        (children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_,
                         false);
 #else
                 children_[i]->applyInjGroupControl(InjectionSpecification::RATE,
-                        (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
+                        (children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_,
                         false);
 #endif
             }
@@ -600,7 +601,7 @@ namespace Opm
             if (phaseUsage().phase_used[BlackoilPhases::Vapour]) {
                 total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Vapour);
             }
-            
+            const double total_reinjected = - total_produced; // Production negative, injection positive
             const double my_guide_rate = injectionGuideRate(true);
             for (size_t i = 0; i < children_.size(); ++i) {
                 // Apply for all children. 
@@ -608,7 +609,7 @@ namespace Opm
                 // as that would check if we're under group control, something we're not.
                 const double children_guide_rate = children_[i]->injectionGuideRate(true);
                 children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
-                        (children_guide_rate / my_guide_rate) * total_produced * injSpec().voidage_replacment_fraction_,
+                        (children_guide_rate / my_guide_rate) * total_reinjected * injSpec().voidage_replacment_fraction_,
                         false);
             }
             
@@ -863,7 +864,7 @@ namespace Opm
             return;
         }
         // We're a producer, so we need to negate the input
-        double ntarget = target;
+        double ntarget = -target;
         
         double distr[3] = { 0.0, 0.0, 0.0 };
         const int* phase_pos = phaseUsage().phase_pos;

opm/core/wells/WellsManager.cpp

@@ -226,6 +226,14 @@ namespace Opm
 
 
 
+    /// Construct from existing wells object.
+    WellsManager::WellsManager(struct Wells* W)
+        : w_(clone_wells(W))
+    {
+    }
+
+
+
     /// Construct wells from deck.
     WellsManager::WellsManager(const Opm::EclipseGridParser& deck,
                                const UnstructuredGrid& grid,

opm/core/wells/WellsManager.hpp

@@ -43,6 +43,13 @@ namespace Opm
 	/// Default constructor -- no wells.
     WellsManager();
 
+    /// Construct from existing wells object.
+    /// WellsManager is not properly initialised in the sense that the logic to
+    /// manage control switching does not exist.
+    ///
+    /// @param[in] W Existing wells object.
+    WellsManager(struct Wells* W);
+
 	/// Construct from input deck and grid.
 	/// The permeability argument may be zero if the input contain
 	/// well productivity indices, otherwise it must be given in

tests/Makefile.am

@@ -1,6 +1,6 @@
 AM_CPPFLAGS =					\
 -I$(top_srcdir)					\
-$(OPM_BOOST_CPPFLAGS)				\
+$(OPM_BOOST_CPPFLAGS) ${SUPERLU_CPPFLAGS}	\
 $(ERT_CPPFLAGS)
 
 AM_LDFLAGS = $(OPM_BOOST_LDFLAGS) $(ERT_LDFLAGS)
@@ -24,6 +24,7 @@ test_read_grid					\
 test_read_vag					\
 test_readpolymer				\
 test_sf2p					\
+test_wells					\
 test_writeVtkData				\
 unit_test
 
@@ -60,6 +61,9 @@ test_sf2p_SOURCES         = test_sf2p.cpp
 
 test_writeVtkData_SOURCES = test_writeVtkData.cpp
 
+test_wells_SOURCES        = test_wells.cpp
+test_wells_LDADD          = $(LDADD) $(BOOST_UNIT_TEST_FRAMEWORK_LIB)
+
 unit_test_SOURCES         = unit_test.cpp
 
 

tests/test_wells.cpp

@@ -0,0 +1,205 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE  // Suppress own messages when throw()ing
+
+#define BOOST_TEST_MODULE WellsModuleTest
+#include <boost/test/unit_test.hpp>
+
+#include <opm/core/newwells.h>
+
+#include <iostream>
+#include <vector>
+#include <boost/shared_ptr.hpp>
+
+BOOST_AUTO_TEST_CASE(Construction)
+{
+    const int nphases = 2;
+    const int nwells  = 2;
+    const int nperfs  = 2;
+
+    boost::shared_ptr<Wells> W(create_wells(nphases, nwells, nperfs),
+                                destroy_wells);
+
+    if (W) {
+        int          cells[] = { 0, 9 };
+        double       WI      = 1.0;
+        const double ifrac[] = { 1.0, 0.0 };
+
+        const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0],
+                                  &WI, "INJECTOR", W.get());
+
+        const double pfrac[] = { 0.0, 0.0 };
+        const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1],
+                                  &WI, "PRODUCER", W.get());
+
+        if (ok0 && ok1) {
+            BOOST_CHECK_EQUAL(W->number_of_phases, nphases);
+            BOOST_CHECK_EQUAL(W->number_of_wells , nwells );
+
+            BOOST_CHECK_EQUAL(W->well_connpos[0], 0);
+            BOOST_CHECK_EQUAL(W->well_connpos[1], 1);
+            BOOST_CHECK_EQUAL(W->well_connpos[W->number_of_wells], nperfs);
+
+            BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[0]], cells[0]);
+            BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[1]], cells[1]);
+
+            BOOST_CHECK_EQUAL(W->WI[W->well_connpos[0]], WI);
+            BOOST_CHECK_EQUAL(W->WI[W->well_connpos[1]], WI);
+
+            using std::string;
+            BOOST_CHECK_EQUAL(string(W->name[0]), string("INJECTOR"));
+            BOOST_CHECK_EQUAL(string(W->name[1]), string("PRODUCER"));
+        }
+    }
+}
+
+
+BOOST_AUTO_TEST_CASE(Controls)
+{
+    const int nphases = 2;
+    const int nwells  = 1;
+    const int nperfs  = 2;
+
+    boost::shared_ptr<Wells> W(create_wells(nphases, nwells, nperfs),
+                               destroy_wells);
+
+    if (W) {
+        int          cells[] = { 0  , 9   };
+        double       WI   [] = { 1.0, 1.0 };
+        const double ifrac[] = { 1.0, 0.0 };
+
+        const bool ok = add_well(INJECTOR, 0.0, nperfs, &ifrac[0], &cells[0],
+                                 &WI[0], "INJECTOR", W.get());
+
+        if (ok) {
+            const double distr[] = { 1.0, 0.0 };
+            const bool   ok1     = append_well_controls(BHP, 1, &distr[0],
+                                                        0, W.get());
+            const bool   ok2     = append_well_controls(SURFACE_RATE, 1,
+                                                        &distr[0], 0, W.get());
+
+            if (ok1 && ok2) {
+                WellControls* ctrls = W->ctrls[0];
+
+                BOOST_CHECK_EQUAL(ctrls->num    ,  2);
+                BOOST_CHECK_EQUAL(ctrls->current, -1);
+
+                set_current_control(0, 0, W.get());
+                BOOST_CHECK_EQUAL(ctrls->current, 0);
+
+                set_current_control(0, 1, W.get());
+                BOOST_CHECK_EQUAL(ctrls->current, 1);
+
+                BOOST_CHECK_EQUAL(ctrls->type[0], BHP);
+                BOOST_CHECK_EQUAL(ctrls->type[1], SURFACE_RATE);
+
+                BOOST_CHECK_EQUAL(ctrls->target[0], 1.0);
+                BOOST_CHECK_EQUAL(ctrls->target[1], 1.0);
+            }
+        }
+    }
+}
+
+
+BOOST_AUTO_TEST_CASE(Copy)
+{
+    const int nphases = 2;
+    const int nwells  = 2;
+    const int nperfs  = 2;
+
+    boost::shared_ptr<Wells> W1(create_wells(nphases, nwells, nperfs),
+                                destroy_wells);
+    boost::shared_ptr<Wells> W2;
+
+    if (W1) {
+        int          cells[] = { 0, 9 };
+        const double WI      = 1.0;
+        const double ifrac[] = { 1.0, 0.0 };
+
+        const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0],
+                                  &WI, "INJECTOR", W1.get());
+
+        const double pfrac[] = { 0.0, 0.0 };
+        const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1],
+                                  &WI, "PRODUCER", W1.get());
+
+        bool ok = ok0 && ok1;
+        for (int w = 0; ok && (w < W1->number_of_wells); ++w) {
+            const double distr[] = { 1.0, 0.0 };
+            const bool   okc1     = append_well_controls(BHP, 1, &distr[0],
+                                                         w, W1.get());
+            const bool   okc2     = append_well_controls(SURFACE_RATE, 1,
+                                                         &distr[0], w,
+                                                         W1.get());
+
+            ok = okc1 && okc2;
+        }
+
+        if (ok) {
+            W2.reset(clone_wells(W1.get()), destroy_wells);
+        }
+    }
+
+    if (W2) {
+        BOOST_CHECK_EQUAL(W2->number_of_phases, W1->number_of_phases);
+        BOOST_CHECK_EQUAL(W2->number_of_wells , W1->number_of_wells );
+        BOOST_CHECK_EQUAL(W2->well_connpos[0] , W1->well_connpos[0] );
+
+        for (int w = 0; w < W1->number_of_wells; ++w) {
+            using std::string;
+            BOOST_CHECK_EQUAL(string(W2->name[w]), string(W1->name[w]));
+            BOOST_CHECK_EQUAL(       W2->type[w] ,        W1->type[w] );
+
+            BOOST_CHECK_EQUAL(W2->well_connpos[w + 1],
+                              W1->well_connpos[w + 1]);
+
+            for (int j = W1->well_connpos[w];
+                j < W1->well_connpos[w + 1]; ++j) {
+                BOOST_CHECK_EQUAL(W2->well_cells[j], W1->well_cells[j]);
+                BOOST_CHECK_EQUAL(W2->WI        [j], W1->WI        [j]);
+            }
+
+            BOOST_CHECK(W1->ctrls[w] != 0);
+            BOOST_CHECK(W2->ctrls[w] != 0);
+
+            WellControls* c1 = W1->ctrls[w];
+            WellControls* c2 = W2->ctrls[w];
+
+            BOOST_CHECK_EQUAL(c2->num    , c1->num    );
+            BOOST_CHECK_EQUAL(c2->current, c1->current);
+
+            for (int c = 0; c < c1->num; ++c) {
+                BOOST_CHECK_EQUAL(c2->type  [c], c1->type  [c]);
+                BOOST_CHECK_EQUAL(c2->target[c], c1->target[c]);
+
+                for (int p = 0; p < W1->number_of_phases; ++p) {
+                    BOOST_CHECK_EQUAL(c2->distr[c*W1->number_of_phases + p],
+                                      c1->distr[c*W1->number_of_phases + p]);
+                }
+            }
+        }
+    }
+}