Rename directory opm/core/tof -> opm/core/flowdiagnostics.

opm/core/flowdiagnostics/AnisotropicEikonal.cpp

@@ -0,0 +1,461 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/core/flowdiagnostics/AnisotropicEikonal.hpp>
+#include <opm/core/grid/GridUtilities.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/utility/RootFinders.hpp>
+
+#if BOOST_HEAP_AVAILABLE
+
+namespace Opm
+{
+
+    namespace
+    {
+        /// Euclidean (isotropic) distance.
+        double distanceIso(const double v1[2],
+                           const double v2[2])
+        {
+            const double d[2] = { v2[0] - v1[0], v2[1] - v1[1] };
+            const double dist = std::sqrt(d[0]*d[0] + d[1]*d[1]);
+            return dist;
+        }
+
+        /// Anisotropic distance with respect to a metric g.
+        /// If d = v2 - v1, the distance is sqrt(d^T g d).
+        double distanceAniso(const double v1[2],
+                             const double v2[2],
+                             const double g[4])
+        {
+            const double d[2] = { v2[0] - v1[0], v2[1] - v1[1] };
+            const double dist = std::sqrt(+ g[0] * d[0] * d[0]
+                                          + g[1] * d[0] * d[1]
+                                          + g[2] * d[1] * d[0]
+                                          + g[3] * d[1] * d[1]);
+            return dist;
+        }
+    } // anonymous namespace
+
+
+
+
+
+
+    /// Construct solver.
+    /// \param[in] grid      A 2d grid.
+    AnisotropicEikonal2d::AnisotropicEikonal2d(const UnstructuredGrid& grid)
+        : grid_(grid),
+          safety_factor_(1.2)
+    {
+        if (grid.dimensions != 2) {
+            OPM_THROW(std::logic_error, "Grid for AnisotropicEikonal2d must be 2d.");
+        }
+        cell_neighbours_ = cellNeighboursAcrossVertices(grid);
+        orderCounterClockwise(grid, cell_neighbours_);
+        computeGridRadius();
+    }
+
+    /// Solve the eikonal equation.
+    /// \param[in]  metric            Array of metric tensors, M, for each cell.
+    /// \param[in]  startcells        Array of cells where u = 0 at the centroid.
+    /// \param[out] solution          Array of solution to the eikonal equation.
+    void AnisotropicEikonal2d::solve(const double* metric,
+                                     const std::vector<int>& startcells,
+                                     std::vector<double>& solution)
+    {
+        // Compute anisotropy ratios to be used by isClose().
+        computeAnisoRatio(metric);
+
+        // The algorithm used is described in J.A. Sethian and A. Vladimirsky,
+        // "Ordered Upwind Methods for Static Hamilton-Jacobi Equations".
+        // Notation in comments is as used in that paper: U is the solution,
+        // and q is the boundary condition. One difference is that we talk about
+        // grid cells instead of mesh points.
+        //
+        // Algorithm summary:
+        // 1. Put all cells in Far. U_i = \inf.
+        // 2. Move the startcells to Accepted. U_i = q(x_i)
+        // 3. Move cells adjacent to startcells to Considered, evaluate
+        //    U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k}
+        // 4. Find the Considered cell with the smallest value: r.
+        // 5. Move cell r to Accepted. Update AcceptedFront.
+        // 6. Recompute the value for all Considered cells within
+        //    distance h * F_2/F1 from x_r. Use min of previous and new.
+        // 7. Move cells adjacent to r from Far to Considered.
+        // 8. If Considered is not empty, go to step 4.
+
+        // 1. Put all cells in Far. U_i = \inf.
+        const int num_cells = grid_.number_of_cells;
+        const double inf = 1e100;
+        solution.clear();
+        solution.resize(num_cells, inf);
+        is_accepted_.clear();
+        is_accepted_.resize(num_cells, false);
+        accepted_front_.clear();
+        considered_.clear();
+        considered_handles_.clear();
+        is_considered_.clear();
+        is_considered_.resize(num_cells, false);
+
+        // 2. Move the startcells to Accepted. U_i = q(x_i)
+        const int num_startcells = startcells.size();
+        for (int ii = 0; ii < num_startcells; ++ii) {
+            is_accepted_[startcells[ii]] = true;
+            solution[startcells[ii]] = 0.0;
+        }
+        accepted_front_.insert(startcells.begin(), startcells.end());
+
+        // 3. Move cells adjacent to startcells to Considered, evaluate
+        //    U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k}
+        for (int ii = 0; ii < num_startcells; ++ii) {
+            const int scell = startcells[ii];
+            const int num_nb = cell_neighbours_[scell].size();
+            for (int nb = 0; nb < num_nb; ++nb) {
+                const int nb_cell = cell_neighbours_[scell][nb];
+                if (!is_accepted_[nb_cell] && !is_considered_[nb_cell]) {
+                    const double value = computeValue(nb_cell, metric, solution.data());
+                    pushConsidered(std::make_pair(value, nb_cell));
+                }
+            }
+        }
+
+        while (!considered_.empty()) {
+            // 4. Find the Considered cell with the smallest value: r.
+            const ValueAndCell r = topConsidered();
+            // std::cout << "Accepting cell " << r.second << std::endl;
+
+            // 5. Move cell r to Accepted. Update AcceptedFront.
+            const int rcell = r.second;
+            is_accepted_[rcell] = true;
+            solution[rcell] = r.first;
+            popConsidered();
+            accepted_front_.insert(rcell);
+            for (auto it = accepted_front_.begin(); it != accepted_front_.end();) {
+                // Note that loop increment happens in the body of this loop.
+                const int cell = *it;
+                bool on_front = false;
+                for (auto it2 = cell_neighbours_[cell].begin(); it2 != cell_neighbours_[cell].end(); ++it2) {
+                    if (!is_accepted_[*it2]) {
+                        on_front = true;
+                        break;
+                    }
+                }
+                if (!on_front) {
+                    accepted_front_.erase(it++);
+                } else {
+                    ++it;
+                }
+            }
+
+            // 6. Recompute the value for all Considered cells within
+            //    distance h * F_2/F1 from x_r. Use min of previous and new.
+            for (auto it = considered_.begin(); it != considered_.end(); ++it) {
+                const int ccell = it->second;
+                if (isClose(rcell, ccell)) {
+                    const double value = computeValueUpdate(ccell, metric, solution.data(), rcell);
+                    if (value < it->first) {
+                        // Update value for considered cell.
+                        // Note that as solution values decrease, their
+                        // goodness w.r.t. the heap comparator increase,
+                        // therefore we may safely call the increase()
+                        // modificator below.
+                        considered_.increase(considered_handles_[ccell], std::make_pair(value, ccell));
+                    }
+                }
+            }
+
+            // 7. Move cells adjacent to r from Far to Considered.
+            for (auto it = cell_neighbours_[rcell].begin(); it != cell_neighbours_[rcell].end(); ++it) {
+                const int nb_cell = *it;
+                if (!is_accepted_[nb_cell] && !is_considered_[nb_cell]) {
+                    assert(solution[nb_cell] == inf);
+                    const double value = computeValue(nb_cell, metric, solution.data());
+                    pushConsidered(std::make_pair(value, nb_cell));
+                }
+            }
+
+            // 8. If Considered is not empty, go to step 4.
+        }
+
+    }
+
+
+
+
+
+    bool AnisotropicEikonal2d::isClose(const int c1,
+                                       const int c2) const
+    {
+        const double* v[] = { grid_.cell_centroids + 2*c1,
+                              grid_.cell_centroids + 2*c2 };
+        return distanceIso(v[0], v[1]) < safety_factor_ * aniso_ratio_[c1] * grid_radius_[c1];
+    }
+
+
+
+
+
+    double AnisotropicEikonal2d::computeValue(const int cell,
+                                              const double* metric,
+                                              const double* solution) const
+    {
+        // std::cout << "++++ computeValue(), cell = " << cell << std::endl;
+        const auto& nbs = cell_neighbours_[cell];
+        const int num_nbs = nbs.size();
+        const double inf = 1e100;
+        double val = inf;
+        for (int ii = 0; ii < num_nbs; ++ii) {
+            const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] };
+            if (accepted_front_.count(n[0]) && accepted_front_.count(n[1])) {
+                const double cand_val = computeFromTri(cell, n[0], n[1], metric, solution);
+                val = std::min(val, cand_val);
+            }
+        }
+        if (val == inf) {
+            // Failed to find two accepted front nodes adjacent to this,
+            // so we go for a single-neighbour update.
+            for (int ii = 0; ii < num_nbs; ++ii) {
+                if (accepted_front_.count(nbs[ii])) {
+                    const double cand_val = computeFromLine(cell, nbs[ii], metric, solution);
+                    val = std::min(val, cand_val);
+                }
+            }
+        }
+        assert(val != inf);
+        // std::cout << "---> " << val << std::endl;
+        return val;
+    }
+
+
+
+
+
+    double AnisotropicEikonal2d::computeValueUpdate(const int cell,
+                                                    const double* metric,
+                                                    const double* solution,
+                                                    const int new_cell) const
+    {
+        // std::cout << "++++ computeValueUpdate(), cell = " << cell << std::endl;
+        const auto& nbs = cell_neighbours_[cell];
+        const int num_nbs = nbs.size();
+        const double inf = 1e100;
+        double val = inf;
+        for (int ii = 0; ii < num_nbs; ++ii) {
+            const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] };
+            if ((n[0] == new_cell || n[1] == new_cell)
+                && accepted_front_.count(n[0]) && accepted_front_.count(n[1])) {
+                const double cand_val = computeFromTri(cell, n[0], n[1], metric, solution);
+                val = std::min(val, cand_val);
+            }
+        }
+        if (val == inf) {
+            // Failed to find two accepted front nodes adjacent to this,
+            // so we go for a single-neighbour update.
+            for (int ii = 0; ii < num_nbs; ++ii) {
+                if (nbs[ii] == new_cell && accepted_front_.count(nbs[ii])) {
+                    const double cand_val = computeFromLine(cell, nbs[ii], metric, solution);
+                    val = std::min(val, cand_val);
+                }
+            }
+        }
+        // std::cout << "---> " << val << std::endl;
+        return val;
+    }
+
+
+
+
+
+    double AnisotropicEikonal2d::computeFromLine(const int cell,
+                                                 const int from,
+                                                 const double* metric,
+                                                 const double* solution) const
+    {
+        assert(!is_accepted_[cell]);
+        assert(is_accepted_[from]);
+        // Applying the first fundamental form to compute geodesic distance.
+        // Using the metric of 'cell', not 'from'.
+        const double dist = distanceAniso(grid_.cell_centroids + 2 * cell,
+                                          grid_.cell_centroids + 2 * from,
+                                          metric + 4 * cell);
+        return solution[from] + dist;
+    }
+
+
+
+
+
+    struct DistanceDerivative
+    {
+        const double* x1;
+        const double* x2;
+        const double* x;
+        double u1;
+        double u2;
+        const double* g;
+        double operator()(const double theta) const
+        {
+            const double xt[2] = { (1-theta)*x1[0] + theta*x2[0], (1-theta)*x1[1] + theta*x2[1] };
+            const double a[2] = { x[0] - xt[0], x[1] - xt[1] };
+            const double b[2] = { x1[0] - x2[0], x1[1] - x2[1] };
+            const double dQdtheta = 2*(a[0]*b[0]*g[0] + a[0]*b[1]*g[1] + a[1]*b[0]*g[2] + a[1]*b[1]*g[3]);
+            const double val =  u2 - u1 + dQdtheta/(2*distanceAniso(x, xt, g));
+            // std::cout << theta << "   " << val << std::endl;
+            return val;
+        }
+    };
+
+
+
+
+
+    double AnisotropicEikonal2d::computeFromTri(const int cell,
+                                                const int n0,
+                                                const int n1,
+                                                const double* metric,
+                                                const double* solution) const
+    {
+        // std::cout << "====  cell = " << cell << "   n0 = " << n0 << "   n1 = " << n1 << std::endl;
+        assert(!is_accepted_[cell]);
+        assert(is_accepted_[n0]);
+        assert(is_accepted_[n1]);
+        DistanceDerivative dd;
+        dd.x1 = grid_.cell_centroids + 2 * n0;
+        dd.x2 = grid_.cell_centroids + 2 * n1;
+        dd.x = grid_.cell_centroids + 2 * cell;
+        dd.u1 = solution[n0];
+        dd.u2 = solution[n1];
+        dd.g = metric + 4 * cell;
+        int iter = 0;
+        const double theta = RegulaFalsi<ContinueOnError>::solve(dd, 0.0, 1.0, 15, 1e-8, iter);
+        const double xt[2] = { (1-theta)*dd.x1[0] + theta*dd.x2[0],
+                               (1-theta)*dd.x1[1] + theta*dd.x2[1] };
+        const double d1 = distanceAniso(dd.x1, dd.x, dd.g) + solution[n0];
+        const double d2 = distanceAniso(dd.x2, dd.x, dd.g) + solution[n1];
+        const double dt = distanceAniso(xt, dd.x, dd.g) + (1-theta)*solution[n0] + theta*solution[n1];
+        return std::min(d1, std::min(d2, dt));
+    }
+
+
+
+
+
+    const AnisotropicEikonal2d::ValueAndCell& AnisotropicEikonal2d::topConsidered() const
+    {
+        return considered_.top();
+    }
+
+
+
+
+
+    void AnisotropicEikonal2d::pushConsidered(const ValueAndCell& vc)
+    {
+        HeapHandle h = considered_.push(vc);
+        considered_handles_[vc.second] = h;
+        is_considered_[vc.second] = true;
+    }
+
+
+
+
+
+    void AnisotropicEikonal2d::popConsidered()
+    {
+        is_considered_[considered_.top().second] = false;
+        considered_handles_.erase(considered_.top().second);
+        considered_.pop();
+    }
+
+
+
+
+    void AnisotropicEikonal2d::computeGridRadius()
+    {
+        const int num_cells = cell_neighbours_.size();
+        grid_radius_.resize(num_cells);
+        for (int cell = 0; cell < num_cells; ++cell) {
+            double radius = 0.0;
+            const double* v1 = grid_.cell_centroids + 2*cell;
+            const auto& nb = cell_neighbours_[cell];
+            for (auto it = nb.begin(); it != nb.end(); ++it) {
+                const double* v2 = grid_.cell_centroids + 2*(*it);
+                radius = std::max(radius, distanceIso(v1, v2));
+            }
+            grid_radius_[cell] = radius;
+        }
+    }
+
+
+
+
+    void AnisotropicEikonal2d::computeAnisoRatio(const double* metric)
+    {
+        const int num_cells = cell_neighbours_.size();
+        aniso_ratio_.resize(num_cells);
+        for (int cell = 0; cell < num_cells; ++cell) {
+            const double* m = metric + 4*cell;
+            // Find the two eigenvalues from trace and determinant.
+            const double t = m[0] + m[3];
+            const double d = m[0]*m[3] - m[1]*m[2];
+            const double sd = std::sqrt(t*t/4.0 - d);
+            const double eig[2] = { t/2.0 - sd, t/2.0 + sd };
+            // Anisotropy ratio is the max ratio of the eigenvalues.
+            aniso_ratio_[cell] = std::max(eig[0]/eig[1], eig[1]/eig[0]);
+        }
+    }
+
+
+
+
+
+} // namespace Opm
+
+
+#else // BOOST_HEAP_AVAILABLE is false
+
+namespace {
+    const char* AnisotropicEikonal2derrmsg =
+        "\n********************************************************************************\n"
+        "This library has not been compiled with support for the AnisotropicEikonal2d\n"
+        "class, due to too old version of the boost libraries (Boost.Heap from boost\n"
+        "version 1.49 or newer is required.\n"
+        "To use this class you must recompile opm-core on a system with sufficiently new\n"
+        "version of the boost libraries."
+        "\n********************************************************************************\n";
+}
+
+namespace Opm
+{
+
+    AnisotropicEikonal2d::AnisotropicEikonal2d(const UnstructuredGrid&)
+    {
+        OPM_THROW(std::logic_error, AnisotropicEikonal2derrmsg);
+    }
+
+    void AnisotropicEikonal2d::solve(const double*,
+                                     const std::vector<int>&,
+                                     std::vector<double>&)
+    {
+        OPM_THROW(std::logic_error, AnisotropicEikonal2derrmsg);
+    }
+}
+
+#endif // BOOST_HEAP_AVAILABLE

opm/core/flowdiagnostics/AnisotropicEikonal.hpp

@@ -0,0 +1,101 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_ANISOTROPICEIKONAL_HEADER_INCLUDED
+#define OPM_ANISOTROPICEIKONAL_HEADER_INCLUDED
+
+#include <opm/core/utility/SparseTable.hpp>
+#include <vector>
+#include <set>
+#include <map>
+
+#include <boost/version.hpp>
+
+#define BOOST_HEAP_AVAILABLE ((BOOST_VERSION / 100 % 1000) >= 49)
+
+#if BOOST_HEAP_AVAILABLE
+#include <boost/heap/fibonacci_heap.hpp>
+#endif
+
+struct UnstructuredGrid;
+
+namespace Opm
+{
+    /// A solver for the anisotropic eikonal equation:
+    ///    \f[ || \nabla u^T M^{-1}(x) \nabla u || = 1 \qquad x \in \Omega \f]
+    /// where M(x) is a symmetric positive definite matrix.
+    /// The boundary conditions are assumed to be
+    ///    \f[ u(x) = 0 \qquad x \in \partial\Omega \f].
+    class AnisotropicEikonal2d
+    {
+    public:
+        /// Construct solver.
+        /// \param[in] grid      A 2d grid.
+        explicit AnisotropicEikonal2d(const UnstructuredGrid& grid);
+
+        /// Solve the eikonal equation.
+        /// \param[in]  metric            Array of metric tensors, M, for each cell.
+        /// \param[in]  startcells        Array of cells where u = 0 at the centroid.
+        /// \param[out] solution          Array of solution to the eikonal equation.
+        void solve(const double* metric,
+                   const std::vector<int>& startcells,
+                   std::vector<double>& solution);
+    private:
+#if BOOST_HEAP_AVAILABLE
+        // Grid and topology.
+        const UnstructuredGrid& grid_;
+        SparseTable<int> cell_neighbours_;
+
+        // Keep track of accepted cells.
+        std::vector<char> is_accepted_;
+        std::set<int> accepted_front_;
+
+        // Quantities relating to anisotropy.
+        std::vector<double> grid_radius_;
+        std::vector<double> aniso_ratio_;
+        const double safety_factor_;
+
+        // Keep track of considered cells.
+        typedef std::pair<double, int> ValueAndCell;
+        typedef boost::heap::compare<std::greater<ValueAndCell>> Comparator;
+        typedef boost::heap::fibonacci_heap<ValueAndCell, Comparator> Heap;
+        Heap considered_;
+        typedef Heap::handle_type HeapHandle;
+        std::map<int, HeapHandle> considered_handles_;
+        std::vector<char> is_considered_;
+
+        bool isClose(const int c1, const int c2) const;
+        double computeValue(const int cell, const double* metric, const double* solution) const;
+        double computeValueUpdate(const int cell, const double* metric, const double* solution, const int new_cell) const;
+        double computeFromLine(const int cell, const int from, const double* metric, const double* solution) const;
+        double computeFromTri(const int cell, const int n0, const int n1, const double* metric, const double* solution) const;
+
+        const ValueAndCell& topConsidered() const;
+        void pushConsidered(const ValueAndCell& vc);
+        void popConsidered();
+
+        void computeGridRadius();
+        void computeAnisoRatio(const double* metric);
+#endif // BOOST_HEAP_AVAILABLE
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_ANISOTROPICEIKONAL_HEADER_INCLUDED

opm/core/flowdiagnostics/DGBasis.cpp

@@ -0,0 +1,342 @@
+/*
+  Copyright 2013 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"
+#include <opm/core/flowdiagnostics/DGBasis.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <numeric>
+
+namespace Opm
+{
+
+    // ----------------  Methods for class DGBasisInterface ----------------
+
+
+    /// Virtual destructor.
+    DGBasisInterface::~DGBasisInterface()
+    {
+    }
+
+    /// Evaluate function f = sum_i c_i b_i at the point x.
+    /// Note that this function is not virtual, but implemented in
+    /// terms of the virtual functions of the class.
+    /// \param[in] cell          Cell index
+    /// \param[in] coefficients  Coefficients {c_i} for a single cell.
+    /// \param[in] x             Point at which to compute f(x).
+    double DGBasisInterface::evalFunc(const int cell,
+                                      const double* coefficients,
+                                      const double* x) const
+    {
+        bvals_.resize(numBasisFunc());
+        eval(cell, x, &bvals_[0]);
+        return std::inner_product(bvals_.begin(), bvals_.end(), coefficients, 0.0);
+    }
+
+
+
+    // ----------------  Methods for class DGBasisBoundedTotalDegree ----------------
+
+
+    /// Constructor.
+    /// \param[in]  grid    grid on which basis is used (cell-wise)
+    /// \param[in]  degree  polynomial degree of basis
+    DGBasisBoundedTotalDegree::DGBasisBoundedTotalDegree(const UnstructuredGrid& grid,
+                                                         const int degree_arg)
+        : grid_(grid),
+          degree_(degree_arg)
+    {
+        if (grid_.dimensions > 3) {
+            OPM_THROW(std::runtime_error, "Grid dimension must be 1, 2 or 3.");
+        }
+        if (degree_ > 1 || degree_ < 0) {
+            OPM_THROW(std::runtime_error, "Degree must be 0 or 1.");
+        }
+    }
+
+    /// Destructor.
+    DGBasisBoundedTotalDegree::~DGBasisBoundedTotalDegree()
+    {
+    }
+
+    /// The number of basis functions per cell.
+    int DGBasisBoundedTotalDegree::numBasisFunc() const
+    {
+        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:
+            OPM_THROW(std::runtime_error, "Dimensions must be 1, 2 or 3.");
+        }
+    }
+
+    /// The number of space dimensions.
+    int DGBasisBoundedTotalDegree::dimensions() const
+    {
+        return grid_.dimensions;
+    }
+
+    /// The polynomial degree of the basis functions.
+    int DGBasisBoundedTotalDegree::degree() const
+    {
+        return degree_;
+    }
+
+    /// Evaluate all basis functions associated with cell at x,
+    /// writing to f_x. The array f_x must have size equal to
+    /// numBasisFunc().
+    void DGBasisBoundedTotalDegree::eval(const int cell,
+                                         const double* x,
+                                         double* f_x) const
+    {
+        const int dim = 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:
+            OPM_THROW(std::runtime_error, "Maximum degree is 1 for now.");
+        }
+    }
+
+
+    /// Evaluate gradients of all basis functions associated with
+    /// cell at x, writing to grad_f_x. The array grad_f_x must
+    /// have size numBasisFunc() * dimensions().  The dimensions()
+    /// components of the first basis function gradient come
+    /// before the components of the second etc.
+    void DGBasisBoundedTotalDegree::evalGrad(const int /*cell*/,
+                                             const double* /*x*/,
+                                             double* grad_f_x) const
+    {
+        const int dim = dimensions();
+        const int num_basis = numBasisFunc();
+        std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0);
+        if (degree_ == 1) {
+            for (int ix = 0; ix < dim; ++ix) {
+                grad_f_x[dim*(ix + 1) + ix] = 1.0;
+            }
+        }
+    }
+
+    /// Modify basis coefficients to add to the function value.
+    /// A function f = sum_i c_i b_i is assumed, and we change
+    /// it to (f + increment) by modifying the c_i. This is done without
+    /// modifying its gradient.
+    /// \param[in]  increment     Add this value to the function.
+    /// \param[out] coefficients  Coefficients {c_i} for a single cell.
+    void DGBasisBoundedTotalDegree::addConstant(const double increment,
+                                                double* coefficients) const
+    {
+        coefficients[0] += increment;
+    }
+
+    /// Modify basis coefficients to change the function's slope.
+    /// A function f = sum_i c_i b_i is assumed, and we change
+    /// it to a function g with the property that grad g = factor * grad f
+    /// by modifying the c_i. This is done without modifying the average,
+    /// i.e. the integrals of g and f over the cell are the same.
+    /// \param[in]  factor        Multiply gradient by this factor.
+    /// \param[out] coefficients  Coefficients {c_i} for a single cell.
+    void DGBasisBoundedTotalDegree::multiplyGradient(const double factor,
+                                                     double* coefficients) const
+    {
+        const int nb = numBasisFunc();
+        for (int ix = 1; ix < nb; ++ix) {
+            coefficients[ix] *= factor;
+        }
+    }
+
+    /// Compute the average of the function f = sum_i c_i b_i.
+    /// \param[in] coefficients  Coefficients {c_i} for a single cell.
+    double DGBasisBoundedTotalDegree::functionAverage(const double* coefficients) const
+    {
+        return coefficients[0];
+    }
+
+
+
+    // ----------------  Methods for class DGBasisMultilin ----------------
+
+
+    /// Constructor.
+    /// \param[in]  grid    grid on which basis is used (cell-wise)
+    /// \param[in]  degree  polynomial degree of basis
+    DGBasisMultilin::DGBasisMultilin(const UnstructuredGrid& grid,
+                                     const int degree_arg)
+        : grid_(grid),
+          degree_(degree_arg)
+    {
+        if (grid_.dimensions > 3) {
+            OPM_THROW(std::runtime_error, "Grid dimension must be 1, 2 or 3.");
+        }
+        if (degree_ > 1 || degree_ < 0) {
+            OPM_THROW(std::runtime_error, "Degree must be 0 or 1.");
+        }
+    }
+
+    /// Destructor.
+    DGBasisMultilin::~DGBasisMultilin()
+    {
+    }
+
+    /// The number of basis functions per cell.
+    int DGBasisMultilin::numBasisFunc() const
+    {
+        switch (dimensions()) {
+        case 1:
+            return degree_ + 1;
+        case 2:
+            return (degree_ + 1)*(degree_ + 1);
+        case 3:
+            return (degree_ + 1)*(degree_ + 1)*(degree_ + 1);
+        default:
+            OPM_THROW(std::runtime_error, "Dimensions must be 1, 2 or 3.");
+        }
+    }
+
+    /// The number of space dimensions.
+    int DGBasisMultilin::dimensions() const
+    {
+        return grid_.dimensions;
+    }
+
+    /// The polynomial degree of the basis functions.
+    int DGBasisMultilin::degree() const
+    {
+        return degree_;
+    }
+
+    /// Evaluate all basis functions associated with cell at x,
+    /// writing to f_x. The array f_x must have size equal to
+    /// numBasisFunc().
+    void DGBasisMultilin::eval(const int cell,
+                               const double* x,
+                               double* f_x) const
+    {
+        const int dim = dimensions();
+        const int num_basis = numBasisFunc();
+        const double* cc = grid_.cell_centroids + dim*cell;
+        switch (degree_) {
+        case 0:
+            f_x[0] = 1;
+            break;
+        case 1:
+            std::fill(f_x, f_x + num_basis, 1.0);
+            for (int dd = 0; dd < dim; ++dd) {
+                const double f[2] = { 0.5 - x[dd] + cc[dd],  0.5 + x[dd] - cc[dd] };
+                const int divi = 1 << (dim - dd - 1); // { 4, 2, 1 } for 3d, for example.
+                for (int ix = 0; ix < num_basis; ++ix) {
+                    f_x[ix] *= f[(ix/divi) % 2];
+                }
+            }
+            break;
+        default:
+            OPM_THROW(std::runtime_error, "Maximum degree is 1 for now.");
+        }
+    }
+
+
+    /// Evaluate gradients of all basis functions associated with
+    /// cell at x, writing to grad_f_x. The array grad_f_x must
+    /// have size numBasisFunc() * dimensions().  The dimensions()
+    /// components of the first basis function gradient come
+    /// before the components of the second etc.
+    void DGBasisMultilin::evalGrad(const int cell,
+                                   const double* x,
+                                   double* grad_f_x) const
+    {
+        const int dim = dimensions();
+        const int num_basis = numBasisFunc();
+            const double* cc = grid_.cell_centroids + dim*cell;
+            switch (degree_) {
+            case 0:
+                std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0);
+                break;
+            case 1:
+                std::fill(grad_f_x, grad_f_x + num_basis*dim, 1.0);
+                for (int dd = 0; dd < dim; ++dd) {
+                    const double f[2] = { 0.5 - x[dd] + cc[dd],  0.5 + x[dd] - cc[dd] };
+                    const double fder[2] = { -1.0,  1.0 };
+                    const int divi = 1 << (dim - dd - 1); // { 4, 2, 1 } for 3d, for example.
+                    for (int ix = 0; ix < num_basis; ++ix) {
+                        const int ind = (ix/divi) % 2;
+                        for (int dder = 0; dder < dim; ++dder) {
+                            grad_f_x[ix*dim + dder] *= (dder == dd ? fder[ind] : f[ind]);
+                        }
+                    }
+                }
+                break;
+            default:
+                OPM_THROW(std::runtime_error, "Maximum degree is 1 for now.");
+            }
+    }
+
+    /// Modify basis coefficients to add to the function value.
+    /// A function f = sum_i c_i b_i is assumed, and we change
+    /// it to (f + increment) by modifying the c_i. This is done without
+    /// modifying its gradient.
+    /// \param[in]  increment     Add this value to the function.
+    /// \param[out] coefficients  Coefficients {c_i} for a single cell.
+    void DGBasisMultilin::addConstant(const double increment,
+                                      double* coefficients) const
+    {
+        const int nb = numBasisFunc();
+        const double term = increment/double(nb);
+        for (int ix = 0; ix < nb; ++ix) {
+            coefficients[ix] += term;
+        }
+    }
+
+    /// Modify basis coefficients to change the function's slope.
+    /// A function f = sum_i c_i b_i is assumed, and we change
+    /// it to a function g with the property that grad g = factor * grad f
+    /// by modifying the c_i. This is done without modifying the average,
+    /// i.e. the integrals of g and f over the cell are the same.
+    /// \param[in]  factor        Multiply gradient by this factor.
+    /// \param[out] coefficients  Coefficients {c_i} for a single cell.
+    void DGBasisMultilin::multiplyGradient(const double factor,
+                                           double* coefficients) const
+    {
+        const int nb = numBasisFunc();
+        const double aver = functionAverage(coefficients);
+        for (int ix = 0; ix < nb; ++ix) {
+            coefficients[ix] = factor*(coefficients[ix] - aver) + aver;
+        }
+    }
+
+    /// Compute the average of the function f = sum_i c_i b_i.
+    /// \param[in] coefficients  Coefficients {c_i} for a single cell.
+    double DGBasisMultilin::functionAverage(const double* coefficients) const
+    {
+        const int nb = numBasisFunc();
+        return std::accumulate(coefficients, coefficients + nb, 0.0)/double(nb);
+    }
+
+} // namespace Opm

opm/core/flowdiagnostics/DGBasis.hpp

@@ -0,0 +1,259 @@
+/*
+  Copyright 2013 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_DGBASIS_HEADER_INCLUDED
+#define OPM_DGBASIS_HEADER_INCLUDED
+
+#include <vector>
+
+struct UnstructuredGrid;
+
+namespace Opm
+{
+
+    /// Base class for Discontinuous Galerkin bases, intended for time-of-flight computations.
+    class DGBasisInterface
+    {
+    public:
+        /// Virtual destructor.
+        virtual ~DGBasisInterface();
+
+        /// The number of basis functions per cell.
+        virtual int numBasisFunc() const = 0;
+
+        /// The number of space dimensions.
+        virtual int dimensions() const = 0;
+
+        /// The polynomial degree of the basis functions.
+        virtual int degree() const = 0;
+
+        /// Evaluate all basis functions associated with cell at x,
+        /// writing to f_x. The array f_x must have size equal to
+        /// numBasisFunc().
+        virtual void eval(const int cell,
+                          const double* x,
+                          double* f_x) const = 0;
+
+        /// Evaluate gradients of all basis functions associated with
+        /// cell at x, writing to grad_f_x. The array grad_f_x must
+        /// have size numBasisFunc() * dimensions().  The dimensions()
+        /// components of the first basis function gradient come
+        /// before the components of the second etc.
+        virtual void evalGrad(const int cell,
+                              const double* x,
+                              double* grad_f_x) const = 0;
+
+        /// Modify basis coefficients to add to the function value.
+        /// A function f = sum_i c_i b_i is assumed, and we change
+        /// it to (f + increment) by modifying the c_i. This is done without
+        /// modifying its gradient.
+        /// \param[in]  increment     Add this value to the function.
+        /// \param[out] coefficients  Coefficients {c_i} for a single cell.
+        virtual void addConstant(const double increment,
+                                 double* coefficients) const = 0;
+
+        /// Modify basis coefficients to change the function's slope.
+        /// A function f = sum_i c_i b_i is assumed, and we change
+        /// it to a function g with the property that grad g = factor * grad f
+        /// by modifying the c_i. This is done without modifying the average,
+        /// i.e. the integrals of g and f over the cell are the same.
+        /// \param[in]  factor        Multiply gradient by this factor.
+        /// \param[out] coefficients  Coefficients {c_i} for a single cell.
+        virtual void multiplyGradient(const double factor,
+                                      double* coefficients) const = 0;
+
+        /// Evaluate function f = sum_i c_i b_i at the point x.
+        /// Note that this function is not virtual, but implemented in
+        /// terms of the virtual functions of the class.
+        /// \param[in] cell          Cell index
+        /// \param[in] coefficients  Coefficients {c_i} for a single cell.
+        /// \param[in] x             Point at which to compute f(x).
+        double evalFunc(const int cell,
+                        const double* coefficients,
+                        const double* x) const;
+
+        /// Compute the average of the function f = sum_i c_i b_i.
+        /// \param[in] coefficients  Coefficients {c_i} for a single cell.
+        virtual double functionAverage(const double* coefficients) const = 0;
+
+    private:
+        mutable std::vector<double> bvals_; // For evalFunc().
+    };
+
+
+
+
+
+    /// A class providing discontinuous Galerkin basis functions
+    /// of bounded total degree.
+    ///
+    /// The basis functions are the following for each cell (example for 3d):
+    ///     Degree 0: 1.
+    ///     Degree 1: 1, x - xc, y - yc, z - zc
+    /// where (xc, yc, zc) are the coordinates of the cell centroid.
+    /// Further degrees await development.
+    class DGBasisBoundedTotalDegree : public DGBasisInterface
+    {
+    public:
+        /// Constructor.
+        /// \param[in]  grid    grid on which basis is used (cell-wise)
+        /// \param[in]  degree  polynomial degree of basis
+        DGBasisBoundedTotalDegree(const UnstructuredGrid& grid, const int degree);
+
+        /// Destructor.
+        virtual ~DGBasisBoundedTotalDegree();
+
+        /// The number of basis functions per cell.
+        virtual int numBasisFunc() const;
+
+        /// The number of space dimensions.
+        virtual int dimensions() const;
+
+        /// The polynomial degree of the basis functions.
+        virtual int degree() const;
+
+        /// Evaluate all basis functions associated with cell at x,
+        /// writing to f_x. The array f_x must have size equal to
+        /// numBasisFunc().
+        virtual void eval(const int cell,
+                          const double* x,
+                          double* f_x) const;
+
+        /// Evaluate gradients of all basis functions associated with
+        /// cell at x, writing to grad_f_x. The array grad_f_x must
+        /// have size numBasisFunc() * dimensions().  The dimensions()
+        /// components of the first basis function gradient come
+        /// before the components of the second etc.
+        virtual void evalGrad(const int cell,
+                              const double* x,
+                              double* grad_f_x) const;
+
+        /// Modify basis coefficients to add to the function value.
+        /// A function f = sum_i c_i b_i is assumed, and we change
+        /// it to (f + increment) by modifying the c_i. This is done without
+        /// modifying its gradient.
+        /// \param[in]  increment     Add this value to the function.
+        /// \param[out] coefficients  Coefficients {c_i} for a single cell.
+        virtual void addConstant(const double increment,
+                                 double* coefficients) const;
+
+        /// Modify basis coefficients to change the function's slope.
+        /// A function f = sum_i c_i b_i is assumed, and we change
+        /// it to a function g with the property that grad g = factor * grad f
+        /// by modifying the c_i. This is done without modifying the average,
+        /// i.e. the integrals of g and f over the cell are the same.
+        /// \param[in]  factor        Multiply gradient by this factor.
+        /// \param[out] coefficients  Coefficients {c_i} for a single cell.
+        virtual void multiplyGradient(const double factor,
+                                      double* coefficients) const;
+
+        /// Compute the average of the function f = sum_i c_i b_i.
+        /// \param[in] coefficients  Coefficients {c_i} for a single cell.
+        virtual double functionAverage(const double* coefficients) const;
+
+    private:
+        const UnstructuredGrid& grid_;
+        const int degree_;
+    };
+
+
+
+
+    /// A class providing discontinuous Galerkin basis functions of
+    /// multi-degree 1 (bilinear or trilinear functions).
+    ///
+    /// The basis functions for a cell are the following
+    ///     Degree 0: 1.
+    /// (for 2 dims:)
+    ///         (Bi)degree 1: (x-)(y-), (x-)(y+), (x+)(y-), (x+)(y+)
+    ///         where (x-) = (1/2 - x + xc), (x+) = (1/2 + x - xc)
+    ///         and xc is the x-coordinate of the cell centroid.
+    ///         Similar for (y-), (y+).
+    class DGBasisMultilin : public DGBasisInterface
+    {
+    public:
+        /// Constructor.
+        /// \param[in]  grid    grid on which basis is used (cell-wise)
+        /// \param[in]  degree  polynomial degree of basis (in each coordinate)
+        DGBasisMultilin(const UnstructuredGrid& grid, const int degree);
+
+        /// Destructor.
+        virtual ~DGBasisMultilin();
+
+        /// The number of basis functions per cell.
+        virtual int numBasisFunc() const;
+
+        /// The number of space dimensions.
+        virtual int dimensions() const;
+
+        /// The polynomial degree of the basis functions.
+        virtual int degree() const;
+
+        /// Evaluate all basis functions associated with cell at x,
+        /// writing to f_x. The array f_x must have size equal to
+        /// numBasisFunc().
+        virtual void eval(const int cell,
+                          const double* x,
+                          double* f_x) const;
+
+        /// Evaluate gradients of all basis functions associated with
+        /// cell at x, writing to grad_f_x. The array grad_f_x must
+        /// have size numBasisFunc() * dimensions().  The dimensions()
+        /// components of the first basis function gradient come
+        /// before the components of the second etc.
+        virtual void evalGrad(const int cell,
+                              const double* x,
+                              double* grad_f_x) const;
+
+        /// Modify basis coefficients to add to the function value.
+        /// A function f = sum_i c_i b_i is assumed, and we change
+        /// it to (f + increment) by modifying the c_i. This is done without
+        /// modifying its gradient.
+        /// \param[in]  increment     Add this value to the function.
+        /// \param[out] coefficients  Coefficients {c_i} for a single cell.
+        virtual void addConstant(const double increment,
+                                 double* coefficients) const;
+
+        /// Modify basis coefficients to change the function's slope.
+        /// A function f = sum_i c_i b_i is assumed, and we change
+        /// it to a function g with the property that grad g = factor * grad f
+        /// by modifying the c_i. This is done without modifying the average,
+        /// i.e. the integrals of g and f over the cell are the same.
+        /// \param[in]  factor        Multiply gradient by this factor.
+        /// \param[out] coefficients  Coefficients {c_i} for a single cell.
+        virtual void multiplyGradient(const double factor,
+                                      double* coefficients) const;
+
+        /// Compute the average of the function f = sum_i c_i b_i.
+        /// \param[in] coefficients  Coefficients {c_i} for a single cell.
+        virtual double functionAverage(const double* coefficients) const;
+
+    private:
+        const UnstructuredGrid& grid_;
+        const int degree_;
+
+    };
+
+
+
+
+} // namespace Opm
+
+
+#endif // OPM_DGBASIS_HEADER_INCLUDED

opm/core/flowdiagnostics/FlowDiagnostics.cpp

@@ -0,0 +1,166 @@
+/*
+  Copyright 2015 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/flowdiagnostics/FlowDiagnostics.hpp>
+
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <algorithm>
+#include <numeric>
+
+namespace Opm
+{
+
+
+    /// \brief Compute flow-capacity/storage-capacity based on time-of-flight.
+    ///
+    /// The F-Phi curve is an analogue to the fractional flow curve in a 1D
+    /// displacement. It can be used to compute other interesting diagnostic
+    /// quantities such as the Lorenz coefficient. For a technical description
+    /// see Shavali et al. (SPE 146446), Shook and Mitchell (SPE 124625).
+    ///
+    /// \param[in]  pv    pore volumes of each cell
+    /// \param[in]  ftof  forward (time from injector) time-of-flight values for each cell
+    /// \param[in]  rtof  reverse (time to producer) time-of-flight values for each cell
+    /// \return           a pair of vectors, the first containing F (flow capacity) the second
+    ///                   containing Phi (storage capacity).
+    std::pair<std::vector<double>, std::vector<double>> computeFandPhi(const std::vector<double>& pv,
+                                                                       const std::vector<double>& ftof,
+                                                                       const std::vector<double>& rtof)
+    {
+        if (pv.size() != ftof.size() || pv.size() != rtof.size()) {
+            OPM_THROW(std::runtime_error, "computeFandPhi(): Input vectors must have same size.");
+        }
+
+        // Sort according to total travel time.
+        const int n = pv.size();
+        typedef std::pair<double, double> D2;
+        std::vector<D2> time_and_pv(n);
+        for (int ii = 0; ii < n; ++ii) {
+            time_and_pv[ii].first = ftof[ii] + rtof[ii]; // Total travel time.
+            time_and_pv[ii].second = pv[ii];
+        }
+        std::sort(time_and_pv.begin(), time_and_pv.end());
+
+        // Compute Phi.
+        std::vector<double> Phi(n + 1);
+        Phi[0] = 0.0;
+        for (int ii = 0; ii < n; ++ii) {
+            Phi[ii+1] = time_and_pv[ii].second;
+        }
+        std::partial_sum(Phi.begin(), Phi.end(), Phi.begin());
+        const double vt = Phi.back(); // Total pore volume.
+        for (int ii = 1; ii < n+1; ++ii) { // Note limits of loop.
+            Phi[ii] /= vt; // Normalize Phi.
+        }
+
+        // Compute F.
+        std::vector<double> F(n + 1);
+        F[0] = 0.0;
+        for (int ii = 0; ii < n; ++ii) {
+            F[ii+1] = time_and_pv[ii].second / time_and_pv[ii].first;
+        }
+        std::partial_sum(F.begin(), F.end(), F.begin());
+        const double ft = F.back(); // Total flux.
+        for (int ii = 1; ii < n+1; ++ii) { // Note limits of loop.
+            F[ii] /= ft; // Normalize Phi.
+        }
+
+        return std::make_pair(F, Phi);
+    }
+
+
+
+
+
+    /// \brief Compute the Lorenz coefficient based on the F-Phi curve.
+    ///
+    /// The Lorenz coefficient is a measure of heterogeneity. It is equal
+    /// to twice the area between the F-Phi curve and the F = Phi line.
+    /// The coefficient can vary from zero to one. If the coefficient is
+    /// zero (so the F-Phi curve is a straight line) we have perfect
+    /// piston-like displacement while a coefficient of one indicates
+    /// infinitely heterogenous displacement (essentially no sweep).
+    ///
+    /// Note: The coefficient is analogous to the Gini coefficient of
+    /// economic theory, where the name Lorenz curve is applied to
+    /// what we call the F-Phi curve.
+    ///
+    /// \param[in]  flowcap     flow capacity (F) as from computeFandPhi()
+    /// \param[in]  storagecap  storage capacity (Phi) as from computeFandPhi()
+    /// \return                 the Lorenz coefficient
+    double computeLorenz(const std::vector<double>& flowcap,
+                         const std::vector<double>& storagecap)
+    {
+        if (flowcap.size() != storagecap.size()) {
+            OPM_THROW(std::runtime_error, "computeLorenz(): Input vectors must have same size.");
+        }
+        double integral = 0.0;
+        // Trapezoid quadrature of the curve F(Phi).
+        const int num_intervals = flowcap.size() - 1;
+        for (int ii = 0; ii < num_intervals; ++ii) {
+            const double len = storagecap[ii+1] - storagecap[ii];
+            integral += (flowcap[ii] + flowcap[ii+1]) * len / 2.0;
+        }
+        return 2.0 * (integral - 0.5);
+    }
+
+
+
+
+
+    /// \brief Compute sweep efficiency versus dimensionless time (PVI).
+    ///
+    /// The sweep efficiency is analogue to 1D displacement using the
+    /// F-Phi curve as flux function.
+    ///
+    /// \param[in]  flowcap     flow capacity (F) as from computeFandPhi()
+    /// \param[in]  storagecap  storage capacity (Phi) as from computeFandPhi()
+    /// \return                 a pair of vectors, the first containing Ev (sweep efficiency)
+    ///                         the second containing tD (dimensionless time).
+    std::pair<std::vector<double>, std::vector<double>> computeSweep(const std::vector<double>& flowcap,
+                                                                     const std::vector<double>& storagecap)
+    {
+        if (flowcap.size() != storagecap.size()) {
+            OPM_THROW(std::runtime_error, "computeSweep(): Input vectors must have same size.");
+        }
+
+        // Compute tD and Ev simultaneously,
+        // skipping identical Phi data points.
+        const int n = flowcap.size();
+        std::vector<double> Ev;
+        std::vector<double> tD;
+        tD.reserve(n);
+        Ev.reserve(n);
+        tD.push_back(0.0);
+        Ev.push_back(0.0);
+        for (int ii = 1; ii < n; ++ii) { // Note loop limits.
+            const double fd = flowcap[ii] - flowcap[ii-1];
+            const double sd = storagecap[ii] - storagecap[ii-1];
+            if (fd != 0.0) {
+                tD.push_back(sd/fd);
+                Ev.push_back(storagecap[ii] + (1.0 - flowcap[ii]) * tD.back());
+            }
+        }
+
+        return std::make_pair(Ev, tD);
+    }
+
+
+
+} // namespace Opm

opm/core/flowdiagnostics/FlowDiagnostics.hpp

@@ -0,0 +1,84 @@
+/*
+  Copyright 2015 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_FLOWDIAGNOSTICS_HEADER_INCLUDED
+#define OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED
+
+
+#include <vector>
+#include <utility>
+
+
+namespace Opm
+{
+
+    /// \brief Compute flow-capacity/storage-capacity based on time-of-flight.
+    ///
+    /// The F-Phi curve is an analogue to the fractional flow curve in a 1D
+    /// displacement. It can be used to compute other interesting diagnostic
+    /// quantities such as the Lorenz coefficient. For a technical description
+    /// see Shavali et al. (SPE 146446), Shook and Mitchell (SPE 124625).
+    ///
+    /// \param[in]  pv    pore volumes of each cell
+    /// \param[in]  ftof  forward (time from injector) time-of-flight values for each cell
+    /// \param[in]  rtof  reverse (time to producer) time-of-flight values for each cell
+    /// \return           a pair of vectors, the first containing F (flow capacity) the second
+    ///                   containing Phi (storage capacity).
+    std::pair<std::vector<double>, std::vector<double>>
+    computeFandPhi(const std::vector<double>& pv,
+                   const std::vector<double>& ftof,
+                   const std::vector<double>& rtof);
+
+
+    /// \brief Compute the Lorenz coefficient based on the F-Phi curve.
+    ///
+    /// The Lorenz coefficient is a measure of heterogeneity. It is equal
+    /// to twice the area between the F-Phi curve and the F = Phi line.
+    /// The coefficient can vary from zero to one. If the coefficient is
+    /// zero (so the F-Phi curve is a straight line) we have perfect
+    /// piston-like displacement while a coefficient of one indicates
+    /// infinitely heterogenous displacement (essentially no sweep).
+    ///
+    /// Note: The coefficient is analogous to the Gini coefficient of
+    /// economic theory, where the name Lorenz curve is applied to
+    /// what we call the F-Phi curve.
+    ///
+    /// \param[in]  flowcap     flow capacity (F) as from computeFandPhi()
+    /// \param[in]  storagecap  storage capacity (Phi) as from computeFandPhi()
+    /// \return                 the Lorenz coefficient
+    double computeLorenz(const std::vector<double>& flowcap,
+                         const std::vector<double>& storagecap);
+
+
+    /// \brief Compute sweep efficiency versus dimensionless time (PVI).
+    ///
+    /// The sweep efficiency is analogue to 1D displacement using the
+    /// F-Phi curve as flux function.
+    ///
+    /// \param[in]  flowcap     flow capacity (F) as from computeFandPhi()
+    /// \param[in]  storagecap  storage capacity (Phi) as from computeFandPhi()
+    /// \return                 a pair of vectors, the first containing Ev (sweep efficiency)
+    ///                         the second containing tD (dimensionless time).
+    std::pair<std::vector<double>, std::vector<double>>
+    computeSweep(const std::vector<double>& flowcap,
+                 const std::vector<double>& storagecap);
+
+} // namespace Opm
+
+#endif // OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED

opm/core/flowdiagnostics/TofDiscGalReorder.cpp

@@ -0,0 +1,805 @@
+/*
+  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"
+#include <opm/core/grid/CellQuadrature.hpp>
+#include <opm/core/grid/FaceQuadrature.hpp>
+#include <opm/core/flowdiagnostics/TofDiscGalReorder.hpp>
+#include <opm/core/flowdiagnostics/DGBasis.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/SparseTable.hpp>
+#include <opm/core/utility/VelocityInterpolation.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/linalg/blas_lapack.h>
+
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+#include <iostream>
+
+namespace Opm
+{
+
+
+    /// Construct solver.
+    TofDiscGalReorder::TofDiscGalReorder(const UnstructuredGrid& grid,
+                                         const parameter::ParameterGroup& param)
+        : grid_(grid),
+          use_cvi_(false),
+          use_limiter_(false),
+          limiter_relative_flux_threshold_(1e-3),
+          limiter_method_(MinUpwindAverage),
+          limiter_usage_(DuringComputations),
+          coord_(grid.dimensions),
+          velocity_(grid.dimensions),
+          gauss_seidel_tol_(1e-3)
+    {
+        const int dg_degree = param.getDefault("dg_degree", 0);
+        const bool use_tensorial_basis = param.getDefault("use_tensorial_basis", false);
+        if (use_tensorial_basis) {
+            basis_func_.reset(new DGBasisMultilin(grid_, dg_degree));
+        } else {
+            basis_func_.reset(new DGBasisBoundedTotalDegree(grid_, dg_degree));
+        }
+
+        tracers_ensure_unity_ = param.getDefault("tracers_ensure_unity", true);
+
+        use_cvi_ = param.getDefault("use_cvi", use_cvi_);
+        use_limiter_ = param.getDefault("use_limiter", use_limiter_);
+        if (use_limiter_) {
+            limiter_relative_flux_threshold_ = param.getDefault("limiter_relative_flux_threshold",
+                                                                limiter_relative_flux_threshold_);
+            const std::string limiter_method_str = param.getDefault<std::string>("limiter_method", "MinUpwindAverage");
+            if (limiter_method_str == "MinUpwindFace") {
+                limiter_method_ = MinUpwindFace;
+            } else if (limiter_method_str == "MinUpwindAverage") {
+                limiter_method_ = MinUpwindAverage;
+            } else {
+                OPM_THROW(std::runtime_error, "Unknown limiter method: " << limiter_method_str);
+            }
+            const std::string limiter_usage_str = param.getDefault<std::string>("limiter_usage", "DuringComputations");
+            if (limiter_usage_str == "DuringComputations") {
+                limiter_usage_ = DuringComputations;
+            } else if (limiter_usage_str == "AsPostProcess") {
+                limiter_usage_ = AsPostProcess;
+            } else if (limiter_usage_str == "AsSimultaneousPostProcess") {
+                limiter_usage_ = AsSimultaneousPostProcess;
+            } else {
+                OPM_THROW(std::runtime_error, "Unknown limiter usage spec: " << limiter_usage_str);
+            }
+        }
+        // A note about the use_cvi_ member variable:
+        // In principle, we should not need it, since the choice of velocity
+        // interpolation is made below, but we may need to use higher order
+        // quadrature to exploit CVI, so we store the choice.
+        // An alternative would be to add a virtual method isConstant() to
+        // the VelocityInterpolationInterface.
+        if (use_cvi_) {
+            velocity_interpolation_.reset(new VelocityInterpolationECVI(grid_));
+        } else {
+            velocity_interpolation_.reset(new VelocityInterpolationConstant(grid_));
+        }
+    }
+
+
+
+
+    /// Solve for time-of-flight.
+    void TofDiscGalReorder::solveTof(const double* darcyflux,
+                                     const double* porevolume,
+                                     const double* source,
+                                     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) {
+            // OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src);
+            OPM_MESSAGE("Warning: sources do not sum to zero: " << cum_src);
+        }
+#endif
+        const int num_basis = basis_func_->numBasisFunc();
+        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);
+        orig_jac_.resize(num_basis*num_basis);
+        basis_.resize(num_basis);
+        basis_nb_.resize(num_basis);
+        grad_basis_.resize(num_basis*grid_.dimensions);
+        velocity_interpolation_->setupFluxes(darcyflux);
+        num_tracers_ = 0;
+        num_multicell_ = 0;
+        max_size_multicell_ = 0;
+        max_iter_multicell_ = 0;
+        num_singlesolves_ = 0;
+        reorderAndTransport(grid_, darcyflux);
+        switch (limiter_usage_) {
+        case AsPostProcess:
+            applyLimiterAsPostProcess();
+            break;
+        case AsSimultaneousPostProcess:
+            applyLimiterAsSimultaneousPostProcess();
+            break;
+        case DuringComputations:
+            // Do nothing.
+            break;
+        default:
+            OPM_THROW(std::runtime_error, "Unknown limiter usage choice: " << limiter_usage_);
+        }
+        if (num_multicell_ > 0) {
+            std::cout << num_multicell_ << " multicell blocks with max size "
+                      << max_size_multicell_ << " cells in upto "
+                      << max_iter_multicell_ << " iterations." << std::endl;
+            std::cout << "Average solves per cell (for all cells) was "
+                      << double(num_singlesolves_)/double(grid_.number_of_cells) << std::endl;
+        }
+    }
+
+
+
+
+    /// Solve for time-of-flight and a number of tracers.
+    /// \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]  tracerheads       Table containing one row per tracer, and each
+    ///                               row contains the source cells for that tracer.
+    /// \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.
+    /// \param[out] tracer_coeff      Array of tracer solution coefficients. N*K per cell,
+    ///                               where N is equal to tracerheads.size(). All K coefs
+    ///                               for a tracer are consecutive, and all tracers' coefs
+    ///                               for a cell come before those for the next cell.
+    void TofDiscGalReorder::solveTofTracer(const double* darcyflux,
+                                           const double* porevolume,
+                                           const double* source,
+                                           const SparseTable<int>& tracerheads,
+                                           std::vector<double>& tof_coeff,
+                                           std::vector<double>& tracer_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) {
+            // OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src);
+            OPM_MESSAGE("Warning: sources do not sum to zero: " << cum_src);
+        }
+#endif
+        const int num_basis = basis_func_->numBasisFunc();
+        num_tracers_ = tracerheads.size();
+        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*(num_tracers_ + 1));
+        jac_.resize(num_basis*num_basis);
+        orig_jac_.resize(num_basis*num_basis);
+        basis_.resize(num_basis);
+        basis_nb_.resize(num_basis);
+        grad_basis_.resize(num_basis*grid_.dimensions);
+        velocity_interpolation_->setupFluxes(darcyflux);
+
+        // Set up tracer
+        tracer_coeff.resize(grid_.number_of_cells*num_tracers_*num_basis);
+        std::fill(tracer_coeff.begin(), tracer_coeff.end(), 0.0);
+        if (num_tracers_ > 0) {
+            tracerhead_by_cell_.clear();
+            tracerhead_by_cell_.resize(grid_.number_of_cells, NoTracerHead);
+        }
+        for (int tr = 0; tr < num_tracers_; ++tr) {
+            for (int i = 0; i < tracerheads[tr].size(); ++i) {
+                const int cell = tracerheads[tr][i];
+                basis_func_->addConstant(1.0, &tracer_coeff[cell*num_tracers_*num_basis + tr*num_basis]);
+                tracer_coeff[cell*num_tracers_ + tr] = 1.0;
+                tracerhead_by_cell_[cell] = tr;
+            }
+        }
+
+        tracer_coeff_ = &tracer_coeff[0];
+        num_multicell_ = 0;
+        max_size_multicell_ = 0;
+        max_iter_multicell_ = 0;
+        num_singlesolves_ = 0;
+        reorderAndTransport(grid_, darcyflux);
+        switch (limiter_usage_) {
+        case AsPostProcess:
+            applyLimiterAsPostProcess();
+            break;
+        case AsSimultaneousPostProcess:
+            applyLimiterAsSimultaneousPostProcess();
+            break;
+        case DuringComputations:
+            // Do nothing.
+            break;
+        default:
+            OPM_THROW(std::runtime_error, "Unknown limiter usage choice: " << limiter_usage_);
+        }
+        if (num_multicell_ > 0) {
+            std::cout << num_multicell_ << " multicell blocks with max size "
+                      << max_size_multicell_ << " cells in upto "
+                      << max_iter_multicell_ << " iterations." << std::endl;
+            std::cout << "Average solves per cell (for all cells) was "
+                      << double(num_singlesolves_)/double(grid_.number_of_cells) << std::endl;
+        }
+    }
+
+
+
+
+    void TofDiscGalReorder::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.
+        //
+        // For tracers, the equation is the same, except for the last
+        // term being zero (the one with \phi).
+        //
+        // The rhs_ vector contains a (Fortran ordering) matrix of all
+        // right-hand-sides, first for tof and then (optionally) for
+        // all tracers.
+
+        const int num_basis = basis_func_->numBasisFunc();
+        ++num_singlesolves_;
+
+        std::fill(rhs_.begin(), rhs_.end(), 0.0);
+        std::fill(jac_.begin(), jac_.end(), 0.0);
+
+        // Add cell contributions to res_ and jac_.
+        cellContribs(cell);
+
+        // Add face contributions to res_ and jac_.
+        faceContribs(cell);
+
+        // Solve linear equation.
+        solveLinearSystem(cell);
+
+        // The solution ends up in rhs_, so we must copy it.
+        std::copy(rhs_.begin(), rhs_.begin() + num_basis, tof_coeff_ + num_basis*cell);
+        if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
+            std::copy(rhs_.begin() + num_basis, rhs_.end(), tracer_coeff_ + num_tracers_*num_basis*cell);
+        }
+
+        // Apply limiter.
+        if (basis_func_->degree() > 0 && use_limiter_ && limiter_usage_ == DuringComputations) {
+            applyLimiter(cell, tof_coeff_);
+            if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
+                for (int tr = 0; tr < num_tracers_; ++tr) {
+                    applyTracerLimiter(cell, tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis);
+                }
+            }
+        }
+
+        // Ensure that tracer averages sum to 1.
+        if (num_tracers_ && tracers_ensure_unity_ && tracerhead_by_cell_[cell] == NoTracerHead) {
+            std::vector<double> tr_aver(num_tracers_);
+            double tr_sum = 0.0;
+            for (int tr = 0; tr < num_tracers_; ++tr) {
+                const double* local_basis = tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis;
+                tr_aver[tr] = basis_func_->functionAverage(local_basis);
+                tr_sum += tr_aver[tr];
+            }
+            if (tr_sum == 0.0) {
+                std::cout << "Tracer sum is zero in cell " << cell << std::endl;
+            } else {
+                for (int tr = 0; tr < num_tracers_; ++tr) {
+                    const double increment = tr_aver[tr]/tr_sum - tr_aver[tr];
+                    double* local_basis = tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis;
+                    basis_func_->addConstant(increment, local_basis);
+                }
+            }
+        }
+    }
+
+
+
+
+    void TofDiscGalReorder::cellContribs(const int cell)
+    {
+        const int num_basis = basis_func_->numBasisFunc();
+        const int dim = grid_.dimensions;
+
+        // Compute cell residual contribution.
+        {
+            const int deg_needed = basis_func_->degree();
+            CellQuadrature quad(grid_, cell, deg_needed);
+            for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
+                // Integral of: b_i \phi
+                quad.quadPtCoord(quad_pt, &coord_[0]);
+                basis_func_->eval(cell, &coord_[0], &basis_[0]);
+                const double w = quad.quadPtWeight(quad_pt);
+                for (int j = 0; j < num_basis; ++j) {
+                    // Only adding to the tof rhs.
+                    rhs_[j] += w * basis_[j] * porevolume_[cell] / grid_.cell_volumes[cell];
+                }
+            }
+        }
+
+        // Compute cell jacobian contribution. We use Fortran ordering
+        // for jac_, i.e. rows cycling fastest.
+        {
+            // Even with ECVI velocity interpolation, degree of precision 1
+            // is sufficient for optimal convergence order for DG1 when we
+            // use linear (total degree 1) basis functions.
+            // With bi(tri)-linear basis functions, it still seems sufficient
+            // for convergence order 2, but the solution looks much better and
+            // has significantly lower error with degree of precision 2.
+            // For now, we err on the side of caution, and use 2*degree, even
+            // though this is wasteful for the pure linear basis functions.
+            // const int deg_needed = 2*basis_func_->degree() - 1;
+            const int deg_needed = 2*basis_func_->degree();
+            CellQuadrature quad(grid_, cell, deg_needed);
+            for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
+                // b_i (v \cdot \grad b_j)
+                quad.quadPtCoord(quad_pt, &coord_[0]);
+                basis_func_->eval(cell, &coord_[0], &basis_[0]);
+                basis_func_->evalGrad(cell, &coord_[0], &grad_basis_[0]);
+                velocity_interpolation_->interpolate(cell, &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 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*basis_func_->degree());
+            for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
+                quad.quadPtCoord(quad_pt, &coord_[0]);
+                basis_func_->eval(cell, &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];
+                    }
+                }
+            }
+        }
+    }
+
+
+
+
+    void TofDiscGalReorder::faceContribs(const int cell)
+    {
+        const int num_basis = basis_func_->numBasisFunc();
+
+        // Compute upstream residual 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;
+            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 (flux >= 0.0) {
+                // This is an outflow boundary.
+                continue;
+            }
+            if (upstream_cell < 0) {
+                // This is an outer boundary. Assumed tof = 0 on inflow, so no contribution.
+                // For tracers, a cell with inflow should be marked as a tracer head cell,
+                // and not be modified.
+                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)).
+            // Quadrature degree set to 2*D, since u_h^{ext} varies
+            // with degree D, and b_j too. We assume that the normal
+            // velocity is constant (this assumption may have to go
+            // for higher order than DG1).
+            const double normal_velocity = flux / grid_.face_areas[face];
+            const int deg_needed = 2*basis_func_->degree();
+            FaceQuadrature quad(grid_, face, deg_needed);
+            for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
+                quad.quadPtCoord(quad_pt, &coord_[0]);
+                basis_func_->eval(cell, &coord_[0], &basis_[0]);
+                basis_func_->eval(upstream_cell, &coord_[0], &basis_nb_[0]);
+                const double w = quad.quadPtWeight(quad_pt);
+                // Modify tof rhs
+                const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(),
+                                                               tof_coeff_ + num_basis*upstream_cell, 0.0);
+                for (int j = 0; j < num_basis; ++j) {
+                    rhs_[j] -= w * tof_upstream * normal_velocity * basis_[j];
+                }
+                // Modify tracer rhs
+                if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
+                    for (int tr = 0; tr < num_tracers_; ++tr) {
+                        const double* up_tr_co = tracer_coeff_ + num_tracers_*num_basis*upstream_cell + num_basis*tr;
+                        const double tracer_up = std::inner_product(basis_nb_.begin(), basis_nb_.end(), up_tr_co, 0.0);
+                        for (int j = 0; j < num_basis; ++j) {
+                            rhs_[num_basis*(tr + 1) + j] -= w * tracer_up * normal_velocity * basis_[j];
+                        }
+                    }
+                }
+            }
+        }
+
+        // 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*basis_func_->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]);
+                basis_func_->eval(cell, &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];
+                    }
+                }
+            }
+        }
+    }
+
+
+
+    // This function assumes that jac_ and rhs_ contain the
+    // linear system to be solved. They are stored in orig_jac_
+    // and orig_rhs_, then the system is solved via LAPACK,
+    // overwriting the input data (jac_ and rhs_).
+    void TofDiscGalReorder::solveLinearSystem(const int cell)
+    {
+        MAT_SIZE_T n = basis_func_->numBasisFunc();
+        int num_tracer_to_compute = num_tracers_;
+        if (num_tracers_) {
+            if (tracerhead_by_cell_[cell] != NoTracerHead) {
+                num_tracer_to_compute = 0;
+            }
+        }
+        MAT_SIZE_T nrhs = 1 + num_tracer_to_compute;
+        MAT_SIZE_T lda = n;
+        std::vector<MAT_SIZE_T> piv(n);
+        MAT_SIZE_T ldb = n;
+        MAT_SIZE_T info = 0;
+        orig_jac_ = jac_;
+        orig_rhs_ = rhs_;
+        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 << "    " << orig_jac_[row + n*col];
+                }
+                std::cerr << '\n';
+            }
+            std::cerr << "and b = \n";
+            for (int row = 0; row < n; ++row) {
+                std::cerr << "    " << orig_rhs_[row] << '\n';
+            }
+            OPM_THROW(std::runtime_error, "Lapack error: " << info << " encountered in cell " << cell);
+        }
+    }
+
+
+
+
+    void TofDiscGalReorder::solveMultiCell(const int num_cells, const int* cells)
+    {
+        ++num_multicell_;
+        max_size_multicell_ = std::max(max_size_multicell_, num_cells);
+        // std::cout << "Multiblock solve with " << num_cells << " cells." << std::endl;
+
+        // Using a Gauss-Seidel approach.
+        const int nb = basis_func_->numBasisFunc();
+        double max_delta = 1e100;
+        int num_iter = 0;
+        while (max_delta > gauss_seidel_tol_) {
+            max_delta = 0.0;
+            ++num_iter;
+            for (int ci = 0; ci < num_cells; ++ci) {
+                const int cell = cells[ci];
+                const double tof_before = basis_func_->functionAverage(&tof_coeff_[nb*cell]);
+                solveSingleCell(cell);
+                const double tof_after = basis_func_->functionAverage(&tof_coeff_[nb*cell]);
+                max_delta = std::max(max_delta, std::fabs(tof_after - tof_before));
+            }
+            // std::cout << "Max delta = " << max_delta << std::endl;
+        }
+        max_iter_multicell_ = std::max(max_iter_multicell_, num_iter);
+    }
+
+
+
+
+    void TofDiscGalReorder::applyLimiter(const int cell, double* tof)
+    {
+        switch (limiter_method_) {
+        case MinUpwindFace:
+            applyMinUpwindLimiter(cell, true, tof);
+            break;
+        case MinUpwindAverage:
+            applyMinUpwindLimiter(cell, false, tof);
+            break;
+        default:
+            OPM_THROW(std::runtime_error, "Limiter type not implemented: " << limiter_method_);
+        }
+    }
+
+
+
+
+    void TofDiscGalReorder::applyMinUpwindLimiter(const int cell, const bool face_min, double* tof)
+    {
+        if (basis_func_->degree() != 1) {
+            OPM_THROW(std::runtime_error, "This limiter only makes sense for our DG1 implementation.");
+        }
+
+        // Limiter principles:
+        // 1. Let M be either:
+        //      - the minimum TOF value of all upstream faces,
+        //        evaluated in the upstream cells
+        //        (chosen if face_min is true).
+        //    or:
+        //      - the minimum average TOF value of all upstream cells
+        //        (chosen if face_min is false).
+        //    Then the value at all points in this cell shall be at
+        //    least M. Upstream faces whose flux does not exceed the
+        //    relative flux threshold are not considered for this
+        //    minimum.
+        // 2. The TOF shall not be below zero in any point.
+
+        // Find minimum tof on upstream faces/cells and for this cell.
+        const int num_basis = basis_func_->numBasisFunc();
+        double min_upstream_tof = 1e100;
+        double min_here_tof = 1e100;
+        int num_upstream_faces = 0;
+        const double total_flux = totalFlux(cell);
+        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];
+            }
+            const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_);
+            const bool interior = (upstream_cell >= 0);
+
+            // Find minimum tof in this cell and upstream.
+            // The meaning of minimum upstream tof depends on method.
+            min_here_tof = std::min(min_here_tof, minCornerVal(cell, face));
+            if (upstream) {
+                ++num_upstream_faces;
+                double upstream_tof = 0.0;
+                if (interior) {
+                    if (face_min) {
+                        upstream_tof = minCornerVal(upstream_cell, face);
+                    } else {
+                        upstream_tof = basis_func_->functionAverage(tof_coeff_ + num_basis*upstream_cell);
+                    }
+                }
+                min_upstream_tof = std::min(min_upstream_tof, upstream_tof);
+            }
+        }
+
+        // Compute slope multiplier (limiter).
+        if (num_upstream_faces == 0) {
+            min_upstream_tof = 0.0;
+            min_here_tof = 0.0;
+        }
+        if (min_upstream_tof < 0.0) {
+            min_upstream_tof = 0.0;
+        }
+        const double tof_c = basis_func_->functionAverage(tof_coeff_ + num_basis*cell);
+        double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof);
+        if (tof_c < min_upstream_tof) {
+            // Handle by setting a flat solution.
+            // std::cout << "Trouble in cell " << cell << std::endl;
+            limiter = 0.0;
+            basis_func_->addConstant(min_upstream_tof - tof_c, tof + num_basis*cell);
+        }
+        assert(limiter >= 0.0);
+
+        // Actually do the limiting (if applicable).
+        if (limiter < 1.0) {
+            // std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl;
+            basis_func_->multiplyGradient(limiter, tof + num_basis*cell);
+        } else {
+            // std::cout << "Not applying limiter in cell " << cell << "!" << std::endl;
+        }
+    }
+
+
+
+
+
+
+    void TofDiscGalReorder::applyLimiterAsPostProcess()
+    {
+        // Apply the limiter sequentially to all cells.
+        // This means that a cell's limiting behaviour may be affected by
+        // any limiting applied to its upstream cells.
+        const std::vector<int>& seq = ReorderSolverInterface::sequence();
+        const int nc = seq.size();
+        assert(nc == grid_.number_of_cells);
+        for (int i = 0; i < nc; ++i) {
+            const int cell = seq[i];
+            applyLimiter(cell, tof_coeff_);
+        }
+    }
+
+
+
+
+    void TofDiscGalReorder::applyLimiterAsSimultaneousPostProcess()
+    {
+        // Apply the limiter simultaneously to all cells.
+        // This means that each cell is limited independently from all other cells,
+        // we write the resulting dofs to a new array instead of writing to tof_coeff_.
+        // Afterwards we copy the results back to tof_coeff_.
+        const int num_basis = basis_func_->numBasisFunc();
+        std::vector<double> tof_coeffs_new(tof_coeff_, tof_coeff_ + num_basis*grid_.number_of_cells);
+        for (int c = 0; c < grid_.number_of_cells; ++c) {
+            applyLimiter(c, &tof_coeffs_new[0]);
+        }
+        std::copy(tof_coeffs_new.begin(), tof_coeffs_new.end(), tof_coeff_);
+    }
+
+
+
+
+    double TofDiscGalReorder::totalFlux(const int cell) const
+    {
+        // Find total upstream/downstream fluxes.
+        double upstream_flux = 0.0;
+        double downstream_flux = 0.0;
+        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) {
+                upstream_flux += flux;
+            } else {
+                downstream_flux += flux;
+            }
+        }
+        // In the presence of sources, significant fluxes may be missing from the computed fluxes,
+        // setting the total flux to the (positive) maximum avoids this: since source is either
+        // inflow or outflow, not both, either upstream_flux or downstream_flux must be correct.
+        return std::max(-upstream_flux, downstream_flux);
+    }
+
+
+
+
+    double TofDiscGalReorder::minCornerVal(const int cell, const int face) const
+    {
+        // Evaluate the solution in all corners.
+        const int dim = grid_.dimensions;
+        const int num_basis = basis_func_->numBasisFunc();
+        double min_cornerval = 1e100;
+        for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) {
+            const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode];
+            basis_func_->eval(cell, nc, &basis_[0]);
+            const double tof_corner = std::inner_product(basis_.begin(), basis_.end(),
+                                                         tof_coeff_ + num_basis*cell, 0.0);
+            min_cornerval = std::min(min_cornerval, tof_corner);
+        }
+        return min_cornerval;
+    }
+
+
+
+    void TofDiscGalReorder::applyTracerLimiter(const int cell, double* local_coeff)
+    {
+        // Evaluate the solution in all corners of all faces. Extract max and min.
+        const int dim = grid_.dimensions;
+        const int num_basis = basis_func_->numBasisFunc();
+        double min_cornerval = 1e100;
+        double max_cornerval = -1e100;
+        for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
+            const int face = grid_.cell_faces[hface];
+            for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) {
+                const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode];
+                basis_func_->eval(cell, nc, &basis_[0]);
+                const double tracer_corner = std::inner_product(basis_.begin(), basis_.end(),
+                                                                local_coeff, 0.0);
+                min_cornerval = std::min(min_cornerval, tracer_corner);
+                max_cornerval = std::max(min_cornerval, tracer_corner);
+            }
+        }
+        const double average = basis_func_->functionAverage(local_coeff);
+        if (average < 0.0 || average > 1.0) {
+            // Adjust average. Flatten gradient.
+            std::fill(local_coeff, local_coeff + num_basis, 0.0);
+            if (average > 1.0) {
+                basis_func_->addConstant(1.0, local_coeff);
+            }
+        } else {
+            // Possibly adjust gradient.
+            double factor = 1.0;
+            if (min_cornerval < 0.0) {
+                factor = average/(average - min_cornerval);
+            }
+            if (max_cornerval > 1.0) {
+                factor = std::min(factor, (1.0 - average)/(max_cornerval - average));
+            }
+            if (factor != 1.0) {
+                basis_func_->multiplyGradient(factor, local_coeff);
+            }
+        }
+    }
+
+
+
+
+} // namespace Opm

opm/core/flowdiagnostics/TofDiscGalReorder.hpp

@@ -0,0 +1,190 @@
+/*
+  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_TOFDISCGALREORDER_HEADER_INCLUDED
+#define OPM_TOFDISCGALREORDER_HEADER_INCLUDED
+
+#include <opm/core/transport/reorder/ReorderSolverInterface.hpp>
+#include <memory>
+#include <vector>
+#include <map>
+#include <ostream>
+
+struct UnstructuredGrid;
+
+namespace Opm
+{
+
+    class IncompPropertiesInterface;
+    class VelocityInterpolationInterface;
+    class DGBasisInterface;
+    namespace parameter { class ParameterGroup; }
+    template <typename T> class SparseTable;
+
+    /// Implements a discontinuous Galerkin solver for
+    /// (single-phase) time-of-flight using reordering.
+    /// The equation solved is:
+    ///     \f[v \cdot \nabla\tau = \phi\f]
+    /// in which \f$ v \f$ is the fluid velocity, \f$ \tau \f$ is time-of-flight and
+    /// \f$ \phi \f$ is the porosity. This is a boundary value problem, and
+    /// \f$ \tau \f$ 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 TofDiscGalReorder : public ReorderSolverInterface
+    {
+    public:
+        /// Construct solver.
+        /// \param[in] grid      A 2d or 3d grid.
+        /// \param[in] param     Parameters for the solver.
+        ///                      The following parameters are accepted (defaults):\n
+        ///   - \c dg_degree (0)                           -- Polynomial degree of basis functions.
+        ///   - \c use_tensorial_basis (false)             -- Use tensor-product basis, interpreting dg_degree as
+        ///                                                   bi/tri-degree not total degree.
+        ///   - \c use_cvi (false)                         -- Use ECVI velocity interpolation.
+        ///   - \c use_limiter (false)                     -- Use a slope limiter. If true, the next three parameters are used.
+        ///   - \c limiter_relative_flux_threshold (1e-3)  -- Ignore upstream fluxes below this threshold,
+        ///                                                   relative to total cell flux.
+        ///   - \c limiter_method ("MinUpwindFace")        -- Limiter method used. Accepted methods are:
+        ///             - MinUpwindFace              -- Limit cell tof to >= inflow face tofs.
+        ///             - MinUpwindAverage           -- Limit cell tof to >= inflow cell average tofs.
+        ///   - \c limiter_usage ("DuringComputations")    -- Usage pattern for limiter. Accepted choices are:
+        ///             - DuringComputations         -- Apply limiter to cells as they are computed,
+        ///                                             so downstream cells' solutions may be affected
+        ///                                             by limiting in upstream cells.
+        ///             - AsPostProcess              -- Apply in dependency order, but only after
+        ///                                             computing (unlimited) solution.
+        ///             - AsSimultaneousPostProcess  -- Apply to each cell independently, using un-
+        ///                                             limited solution in neighbouring cells.
+        TofDiscGalReorder(const UnstructuredGrid& grid,
+                          const parameter::ParameterGroup& param);
+
+
+        /// 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_coeff         Array of time-of-flight solution coefficients.
+        ///                               The values are ordered by cell, meaning that
+        ///                               the K coefficients corresponding to the first
+        ///                               cell come 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,
+                      std::vector<double>& tof_coeff);
+
+        /// Solve for time-of-flight and a number of tracers.
+        /// \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]  tracerheads       Table containing one row per tracer, and each
+        ///                               row contains the source cells for that tracer.
+        /// \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.
+        /// \param[out] tracer_coeff      Array of tracer solution coefficients. N*K per cell,
+        ///                               where N is equal to tracerheads.size(). All K coefs
+        ///                               for a tracer are consecutive, and all tracers' coefs
+        ///                               for a cell come before those for the next cell.
+        void solveTofTracer(const double* darcyflux,
+                            const double* porevolume,
+                            const double* source,
+                            const SparseTable<int>& tracerheads,
+                            std::vector<double>& tof_coeff,
+                            std::vector<double>& tracer_coeff);
+
+    private:
+        virtual void solveSingleCell(const int cell);
+        virtual void solveMultiCell(const int num_cells, const int* cells);
+
+        void cellContribs(const int cell);
+        void faceContribs(const int cell);
+        void solveLinearSystem(const int cell);
+
+    private:
+        // Disable copying and assignment.
+        TofDiscGalReorder(const TofDiscGalReorder&);
+        TofDiscGalReorder& operator=(const TofDiscGalReorder&);
+
+        // Data members
+        const UnstructuredGrid& grid_;
+        std::shared_ptr<VelocityInterpolationInterface> velocity_interpolation_;
+        bool use_cvi_;
+        bool use_limiter_;
+        double limiter_relative_flux_threshold_;
+        enum LimiterMethod { MinUpwindFace, MinUpwindAverage };
+        LimiterMethod limiter_method_;
+        enum LimiterUsage { DuringComputations, AsPostProcess, AsSimultaneousPostProcess };
+        LimiterUsage limiter_usage_;
+        const double* darcyflux_;   // one flux per grid face
+        const double* porevolume_;  // one volume per cell
+        const double* source_;      // one volumetric source term per cell
+        std::shared_ptr<DGBasisInterface> basis_func_;
+        double* tof_coeff_;
+        // For tracers.
+        double* tracer_coeff_;
+        int num_tracers_;
+        enum { NoTracerHead = -1 };
+        std::vector<int> tracerhead_by_cell_;
+        bool tracers_ensure_unity_;
+        // Used by solveSingleCell().
+        std::vector<double> rhs_;   // single-cell right-hand-sides
+        std::vector<double> jac_;   // single-cell jacobian
+        std::vector<double> orig_rhs_;   // single-cell right-hand-sides (copy)
+        std::vector<double> orig_jac_;   // single-cell jacobian (copy)
+        std::vector<double> coord_;
+        mutable std::vector<double> basis_;
+        mutable std::vector<double> basis_nb_;
+        std::vector<double> grad_basis_;
+        std::vector<double> velocity_;
+        int num_singlesolves_;
+        // Used by solveMultiCell():
+        double gauss_seidel_tol_;
+        int num_multicell_;
+        int max_size_multicell_;
+        int max_iter_multicell_;
+
+        // Private methods
+
+        // Apply some limiter, writing to array tof
+        // (will read data from tof_coeff_, it is ok to call
+        //  with tof_coeff as tof argument.
+        void applyLimiter(const int cell, double* tof);
+        void applyMinUpwindLimiter(const int cell, const bool face_min, double* tof);
+        void applyLimiterAsPostProcess();
+        void applyLimiterAsSimultaneousPostProcess();
+        double totalFlux(const int cell) const;
+        double minCornerVal(const int cell, const int face) const;
+
+        // Apply a simple (restrict to [0,1]) limiter.
+        // Intended for tracers.
+        void applyTracerLimiter(const int cell, double* local_coeff);
+    };
+
+} // namespace Opm
+
+#endif // OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED

opm/core/flowdiagnostics/TofReorder.cpp

@@ -0,0 +1,452 @@
+/*
+  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"
+#include <opm/core/flowdiagnostics/TofReorder.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/SparseTable.hpp>
+
+#include <algorithm>
+#include <numeric>
+#include <cmath>
+#include <iostream>
+
+namespace Opm
+{
+
+
+    /// Construct solver.
+    /// \param[in] grid      A 2d or 3d grid.
+    /// \param[in] use_multidim_upwind  If true, use multidimensional tof upwinding.
+    TofReorder::TofReorder(const UnstructuredGrid& grid,
+                           const bool use_multidim_upwind)
+        : grid_(grid),
+          darcyflux_(0),
+          porevolume_(0),
+          source_(0),
+          tof_(0),
+          gauss_seidel_tol_(1e-3),
+          use_multidim_upwind_(use_multidim_upwind)
+    {
+    }
+
+
+
+
+    /// 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 TofReorder::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) {
+            // OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src);
+            OPM_MESSAGE("Warning: 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];
+        if (use_multidim_upwind_) {
+            face_tof_.resize(grid_.number_of_faces);
+            std::fill(face_tof_.begin(), face_tof_.end(), 0.0);
+            face_part_tof_.resize(grid_.face_nodepos[grid_.number_of_faces]);
+            std::fill(face_part_tof_.begin(), face_part_tof_.end(), 0.0);
+        }
+        compute_tracer_ = false;
+        executeSolve();
+    }
+
+
+
+
+    /// Solve for time-of-flight and a number of tracers.
+    /// \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]  tracerheads       Table containing one row per tracer, and each
+    ///                               row contains the source cells for that tracer.
+    /// \param[out] tof               Array of time-of-flight values (1 per cell).
+    /// \param[out] tracer            Array of tracer values. N per cell, where N is
+    ///                               equalt to tracerheads.size().
+    void TofReorder::solveTofTracer(const double* darcyflux,
+                                    const double* porevolume,
+                                    const double* source,
+                                    const SparseTable<int>& tracerheads,
+                                    std::vector<double>& tof,
+                                    std::vector<double>& tracer)
+    {
+        darcyflux_ = darcyflux;
+        porevolume_ = porevolume;
+        source_ = source;
+        const int num_cells = grid_.number_of_cells;
+#ifndef NDEBUG
+        // Sanity check for sources.
+        const double cum_src = std::accumulate(source, source + num_cells, 0.0);
+        if (std::fabs(cum_src) > *std::max_element(source, source + num_cells)*1e-2) {
+            OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src);
+        }
+#endif
+        tof.resize(num_cells);
+        std::fill(tof.begin(), tof.end(), 0.0);
+        tof_ = &tof[0];
+
+        if (use_multidim_upwind_) {
+            face_tof_.resize(grid_.number_of_faces);
+            std::fill(face_tof_.begin(), face_tof_.end(), 0.0);
+            face_part_tof_.resize(grid_.face_nodepos[grid_.number_of_faces]);
+            std::fill(face_part_tof_.begin(), face_part_tof_.end(), 0.0);
+        }
+
+        // Execute solve for tof
+        compute_tracer_ = false;
+        executeSolve();
+
+        // Find the tracer heads (injectors).
+        const int num_tracers = tracerheads.size();
+        tracer.resize(num_cells*num_tracers);
+        std::fill(tracer.begin(), tracer.end(), 0.0);
+        if (num_tracers > 0) {
+            tracerhead_by_cell_.clear();
+            tracerhead_by_cell_.resize(num_cells, NoTracerHead);
+        }
+        for (int tr = 0; tr < num_tracers; ++tr) {
+            for (int i = 0; i < tracerheads[tr].size(); ++i) {
+                const int cell = tracerheads[tr][i];
+                tracer[num_cells * tr + cell] = 1.0;
+                tracerhead_by_cell_[cell] = tr;
+            }
+        }
+
+        // Execute solve for tracers.
+        std::vector<double> fake_pv(num_cells, 0.0);
+        porevolume_ = fake_pv.data();
+        for (int tr = 0; tr < num_tracers; ++tr) {
+            tof_ = tracer.data() + tr * num_cells;
+            compute_tracer_ = true;
+            executeSolve();
+        }
+
+        // Write output tracer data (transposing the computed data).
+        std::vector<double> computed = tracer;
+        for (int cell = 0; cell < num_cells; ++cell) {
+            for (int tr = 0; tr < num_tracers; ++tr) {
+                tracer[num_tracers * cell + tr] = computed[num_cells * tr + cell];
+            }
+        }
+    }
+
+
+
+
+    void TofReorder::executeSolve()
+    {
+        num_multicell_ = 0;
+        max_size_multicell_ = 0;
+        max_iter_multicell_ = 0;
+        reorderAndTransport(grid_, darcyflux_);
+        if (num_multicell_ > 0) {
+            std::cout << num_multicell_ << " multicell blocks with max size "
+                      << max_size_multicell_ << " cells in upto "
+                      << max_iter_multicell_ << " iterations." << std::endl;
+        }
+    }
+
+
+
+    void TofReorder::solveSingleCell(const int cell)
+    {
+        if (use_multidim_upwind_) {
+            solveSingleCellMultidimUpwind(cell);
+            return;
+        }
+        // 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).
+        if (compute_tracer_ && tracerhead_by_cell_[cell] != NoTracerHead) {
+            // This is a tracer head cell, already has solution.
+            return;
+        }
+        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 (flux < 0.0) {
+                // Using tof == 0 on inflow, so we only add a
+                // nonzero contribution if we are on an internal
+                // face.
+                if (other != -1) {
+                    upwind_term += flux*tof_[other];
+                }
+            } else {
+                downwind_flux += flux;
+            }
+        }
+
+        // Compute tof.
+        tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux;
+    }
+
+
+
+
+    void TofReorder::solveSingleCellMultidimUpwind(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 terms (note sign change resulting from
+        // different sign conventions: pos. source is injection,
+        // pos. flux is outflow).
+        double upwind_term = 0.0;
+        double downwind_term_cell_factor = std::max(-source_[cell], 0.0);
+        double downwind_term_face = 0.0;
+        for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
+            int f = grid_.cell_faces[i];
+            double flux;
+            // Compute cell flux
+            if (cell == grid_.face_cells[2*f]) {
+                flux  = darcyflux_[f];
+            } else {
+                flux  =-darcyflux_[f];
+            }
+            // Add flux to upwind_term or downwind_term_[face|cell_factor].
+            if (flux < 0.0) {
+                upwind_term += flux*face_tof_[f];
+            } else if (flux > 0.0) {
+                double fterm, cterm_factor;
+                multidimUpwindTerms(f, cell, fterm, cterm_factor);
+                downwind_term_face += fterm*flux;
+                downwind_term_cell_factor += cterm_factor*flux;
+            }
+        }
+
+        // Compute tof for cell.
+        if (compute_tracer_ && tracerhead_by_cell_[cell] != NoTracerHead) {
+            // Do nothing to the value in this cell, since we are at a tracer head.
+        } else {
+            tof_[cell] = (porevolume_[cell] - upwind_term - downwind_term_face)/downwind_term_cell_factor;
+        }
+
+        // Compute tof for downwind faces.
+        for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
+            int f = grid_.cell_faces[i];
+            const double outflux_f = (grid_.face_cells[2*f] == cell) ? darcyflux_[f] : -darcyflux_[f];
+            if (outflux_f > 0.0) {
+                double fterm, cterm_factor;
+                multidimUpwindTerms(f, cell, fterm, cterm_factor);
+                face_tof_[f] = fterm + cterm_factor*tof_[cell];
+
+                // Combine locally computed (for each adjacent vertex) terms, with uniform weighting.
+                const int* face_nodes_beg = grid_.face_nodes + grid_.face_nodepos[f];
+                const int* face_nodes_end = grid_.face_nodes + grid_.face_nodepos[f + 1];
+                assert((face_nodes_end - face_nodes_beg) == 2 || grid_.dimensions != 2);
+                for (const int* fn_iter = face_nodes_beg; fn_iter < face_nodes_end; ++fn_iter) {
+                    double loc_face_term = 0.0;
+                    double loc_cell_term_factor = 0.0;
+                    const int node_pos = fn_iter - grid_.face_nodes;
+                    localMultidimUpwindTerms(f, cell, node_pos,
+                                             loc_face_term, loc_cell_term_factor);
+                    face_part_tof_[node_pos] = loc_face_term + loc_cell_term_factor * tof_[cell];
+                }
+            }
+        }
+    }
+
+
+
+
+    void TofReorder::solveMultiCell(const int num_cells, const int* cells)
+    {
+        ++num_multicell_;
+        max_size_multicell_ = std::max(max_size_multicell_, num_cells);
+        // std::cout << "Multiblock solve with " << num_cells << " cells." << std::endl;
+
+        // Using a Gauss-Seidel approach.
+        double max_delta = 1e100;
+        int num_iter = 0;
+        while (max_delta > gauss_seidel_tol_) {
+            max_delta = 0.0;
+            ++num_iter;
+            for (int ci = 0; ci < num_cells; ++ci) {
+                const int cell = cells[ci];
+                const double tof_before = tof_[cell];
+                solveSingleCell(cell);
+                max_delta = std::max(max_delta, std::fabs(tof_[cell] - tof_before));
+            }
+            // std::cout << "Max delta = " << max_delta << std::endl;
+        }
+        max_iter_multicell_ = std::max(max_iter_multicell_, num_iter);
+    }
+
+
+
+
+    // Assumes that face_part_tof_[node_pos] is known for all inflow
+    // faces to 'upwind_cell' sharing vertices with 'face'. The index
+    // 'node_pos' is the same as the one used for the grid face-node
+    // connectivity.
+    // Assumes that darcyflux_[face] is != 0.0.
+    // This function returns factors to compute the tof for 'face':
+    //   tof(face) = face_term + cell_term_factor*tof(upwind_cell).
+    // It is not computed here, since these factors are needed to
+    // compute the tof(upwind_cell) itself.
+    void TofReorder::multidimUpwindTerms(const int face,
+                                         const int upwind_cell,
+                                         double& face_term,
+                                         double& cell_term_factor) const
+    {
+        // Implements multidim upwind inspired by
+        // "Multidimensional upstream weighting for multiphase transport on general grids"
+        // by Keilegavlen, Kozdon, Mallison.
+        // However, that article does not give a 3d extension other than noting that using
+        // multidimensional upwinding in the XY-plane and not in the Z-direction may be
+        // a good idea. We have here attempted some generalization, by treating each face-part
+        // (association of a face and a vertex) as possibly influencing all downwind face-parts
+        // of the neighbouring cell that share the same vertex.
+        // The current implementation aims to reproduce 2d results for extruded 3d grids.
+
+        // Combine locally computed (for each adjacent vertex) terms, with uniform weighting.
+        const int* face_nodes_beg = grid_.face_nodes + grid_.face_nodepos[face];
+        const int* face_nodes_end = grid_.face_nodes + grid_.face_nodepos[face + 1];
+        const int num_terms = face_nodes_end - face_nodes_beg;
+        assert(num_terms == 2 || grid_.dimensions != 2);
+        face_term = 0.0;
+        cell_term_factor = 0.0;
+        for (const int* fn_iter = face_nodes_beg; fn_iter < face_nodes_end; ++fn_iter) {
+            double loc_face_term = 0.0;
+            double loc_cell_term_factor = 0.0;
+            localMultidimUpwindTerms(face, upwind_cell, fn_iter - grid_.face_nodes,
+                                     loc_face_term, loc_cell_term_factor);
+            face_term += loc_face_term;
+            cell_term_factor += loc_cell_term_factor;
+        }
+        face_term /= double(num_terms);
+        cell_term_factor /= double(num_terms);
+
+    }
+
+
+
+
+    namespace {
+        double weightFunc(const double w)
+        {
+            // SPU
+            // return 0.0;
+            // TMU
+            return w > 0.0 ? std::min(w, 1.0) :  0.0;
+            // SMU
+            // return w > 0.0 ? w/(1.0 + w) : 0.0;
+        }
+    }
+
+
+
+
+    void TofReorder::localMultidimUpwindTerms(const int face,
+                                              const int upwind_cell,
+                                              const int node_pos,
+                                              double& face_term,
+                                              double& cell_term_factor) const
+    {
+        // Loop over all faces adjacent to the given cell and the
+        // vertex in position node_pos.
+        // If that part's influx is positive, we store it, and also its associated
+        // node position.
+        std::vector<double> influx;
+        std::vector<int> node_pos_influx;
+        influx.reserve(5);
+        node_pos_influx.reserve(5);
+        const int node = grid_.face_nodes[node_pos];
+        for (int hf = grid_.cell_facepos[upwind_cell]; hf < grid_.cell_facepos[upwind_cell + 1]; ++hf) {
+            const int f = grid_.cell_faces[hf];
+            if (f != face) {
+                // Find out if the face 'f' is adjacent to vertex 'node'.
+                const int* f_nodes_beg = grid_.face_nodes + grid_.face_nodepos[f];
+                const int* f_nodes_end = grid_.face_nodes + grid_.face_nodepos[f + 1];
+                const int* pos = std::find(f_nodes_beg, f_nodes_end, node);
+                const int node_pos2 = pos - grid_.face_nodes;
+                const bool is_adj = (pos != f_nodes_end);
+                if (is_adj) {
+                    const int num_parts = f_nodes_end - f_nodes_beg;
+                    const double influx_sign = (grid_.face_cells[2*f] == upwind_cell) ? -1.0 : 1.0;
+                    const double part_influx = influx_sign * darcyflux_[f] / double(num_parts);
+                    if (part_influx > 0.0) {
+                        influx.push_back(part_influx);
+                        node_pos_influx.push_back(node_pos2);
+                    }
+                }
+            }
+        }
+
+        // Now we may compute the weighting of the upwind terms.
+        const int num_parts = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
+        const double outflux_sign = (grid_.face_cells[2*face] == upwind_cell) ? 1.0 : -1.0;
+        const double part_outflux = outflux_sign * darcyflux_[face] / double(num_parts);
+        const double sum_influx = std::accumulate(influx.begin(), influx.end(), 0.0);
+        const double w_factor = weightFunc(sum_influx / part_outflux);
+        const int num_influx = influx.size();
+        std::vector<double> w(num_influx);
+        face_term = 0.0;
+        for (int ii = 0; ii < num_influx; ++ii) {
+            w[ii] = (influx[ii] / sum_influx) * w_factor;
+            face_term += w[ii] * face_part_tof_[node_pos_influx[ii]];
+        }
+        const double sum_w = std::accumulate(w.begin(), w.end(), 0.0);
+        cell_term_factor = 1.0 - sum_w;
+        const double tol = 1e-5;
+        if (cell_term_factor < -tol && cell_term_factor > 1.0 + tol) {
+            OPM_THROW(std::logic_error, "cell_term_factor outside [0,1]: " << cell_term_factor);
+        }
+        cell_term_factor = std::min(std::max(cell_term_factor, 0.0), 1.0);
+        assert(cell_term_factor >= 0.0);
+        assert(cell_term_factor <= 1.0);
+    }
+
+} // namespace Opm

opm/core/flowdiagnostics/TofReorder.hpp

@@ -0,0 +1,119 @@
+/*
+  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_TOFREORDER_HEADER_INCLUDED
+#define OPM_TOFREORDER_HEADER_INCLUDED
+
+#include <opm/core/transport/reorder/ReorderSolverInterface.hpp>
+#include <vector>
+#include <map>
+#include <ostream>
+
+struct UnstructuredGrid;
+
+namespace Opm
+{
+
+    class IncompPropertiesInterface;
+    template <typename T> class SparseTable;
+
+    /// Implements a first-order finite volume solver for
+    /// (single-phase) time-of-flight using reordering.
+    /// The equation solved is:
+    ///     \f[v \cdot \nabla\tau = \phi\f]
+    /// in which \f$ v \f$ is the fluid velocity, \f$ \tau \f$ is time-of-flight and
+    /// \f$ \phi \f$ is the porosity. This is a boundary value problem, and
+    /// \f$ \tau \f$ is specified to be zero on all inflow boundaries.
+    class TofReorder : public ReorderSolverInterface
+    {
+    public:
+        /// Construct solver.
+        /// \param[in] grid      A 2d or 3d grid.
+        /// \param[in] use_multidim_upwind  If true, use multidimensional tof upwinding.
+        TofReorder(const UnstructuredGrid& grid,
+                   const bool use_multidim_upwind = false);
+
+        /// 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);
+
+        /// Solve for time-of-flight and a number of tracers.
+        /// \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]  tracerheads       Table containing one row per tracer, and each
+        ///                               row contains the source cells for that tracer.
+        /// \param[out] tof               Array of time-of-flight values (1 per cell).
+        /// \param[out] tracer            Array of tracer values. N per cell, where N is
+        ///                               equalt to tracerheads.size().
+        void solveTofTracer(const double* darcyflux,
+                            const double* porevolume,
+                            const double* source,
+                            const SparseTable<int>& tracerheads,
+                            std::vector<double>& tof,
+                            std::vector<double>& tracer);
+
+    private:
+        void executeSolve();
+        virtual void solveSingleCell(const int cell);
+        void solveSingleCellMultidimUpwind(const int cell);
+        void assembleSingleCell(const int cell,
+                                std::vector<int>& local_column,
+                                std::vector<double>& local_coefficient,
+                                double& rhs);
+        virtual void solveMultiCell(const int num_cells, const int* cells);
+
+        void multidimUpwindTerms(const int face, const int upwind_cell,
+                                 double& face_term, double& cell_term_factor) const;
+        void localMultidimUpwindTerms(const int face, const int upwind_cell, const int node_pos,
+                                      double& face_term, double& cell_term_factor) const;
+
+    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_;
+        bool compute_tracer_;
+        enum { NoTracerHead = -1 };
+        std::vector<int> tracerhead_by_cell_;
+        // For solveMultiCell():
+        double gauss_seidel_tol_;
+        int num_multicell_;
+        int max_size_multicell_;
+        int max_iter_multicell_;
+        // For multidim upwinding:
+        bool use_multidim_upwind_;
+        std::vector<double> face_tof_;       // For multidim upwind face tofs.
+        std::vector<double> face_part_tof_;  // For multidim upwind face tofs.
+    };
+
+} // namespace Opm
+
+#endif // OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED