Opdated opm-flowdiagnostics to sha: d409aedbf96c996c3886e3fd91abd4d8c3b0c601

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/CMakeLists_files.cmake

@@ -25,9 +25,9 @@
 
 list (APPEND MAIN_SOURCE_FILES
         opm/flowdiagnostics/CellSet.cpp
-        opm/flowdiagnostics/CellSetValues.cpp
         opm/flowdiagnostics/ConnectionValues.cpp
         opm/flowdiagnostics/ConnectivityGraph.cpp
+        opm/flowdiagnostics/DerivedQuantities.cpp
         opm/flowdiagnostics/Solution.cpp
         opm/flowdiagnostics/Toolbox.cpp
         opm/flowdiagnostics/TracerTofSolver.cpp
@@ -39,9 +39,9 @@ list (APPEND MAIN_SOURCE_FILES
 list (APPEND TEST_SOURCE_FILES
         tests/test_assembledconnections.cpp
         tests/test_cellset.cpp
-        tests/test_cellsetvalues.cpp
         tests/test_connectionvalues.cpp
         tests/test_connectivitygraph.cpp
+        tests/test_derivedquantities.cpp
         tests/test_flowdiagnosticstool.cpp
         tests/test_tarjan.cpp
         )
@@ -51,6 +51,7 @@ list (APPEND PUBLIC_HEADER_FILES
         opm/flowdiagnostics/CellSetValues.hpp
         opm/flowdiagnostics/ConnectionValues.hpp
         opm/flowdiagnostics/ConnectivityGraph.hpp
+        opm/flowdiagnostics/DerivedQuantities.hpp
         opm/flowdiagnostics/Solution.hpp
         opm/flowdiagnostics/Toolbox.hpp
         opm/flowdiagnostics/TracerTofSolver.hpp

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/CellSet.cpp

@@ -46,12 +46,22 @@ CellSetID::to_string() const
     return id_;
 }
 
+bool
+CellSetID::operator<(const CellSetID& other) const
+{
+    return id_ < other.id_;
+}
+
 // =====================================================================
 
-void
+CellSet::CellSet(CellSetID id)
-CellSet::identify(CellSetID id)
+    : id_(std::move(id))
+{
+}
+
+CellSet::CellSet(CellSetID id, const std::vector<int>& cells)
+    : id_(std::move(id)), iset_(cells.begin(), cells.end())
 {
-    id_ = std::move(id);
 }
 
 const CellSetID&

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/CellSet.hpp

@@ -23,6 +23,7 @@
 
 #include <string>
 #include <unordered_set>
+#include <vector>
 
 namespace Opm
 {
@@ -41,24 +42,32 @@ namespace FlowDiagnostics
 
         std::string to_string() const;
 
+        bool operator<(const CellSetID& other) const;
+
     private:
         Repr id_;
     };
 
+
+
+
     class CellSet
     {
     private:
         using IndexSet = std::unordered_set<int>;
 
     public:
-        using const_iterator = IndexSet::const_iterator;
+        /// Contruct empty cell set, use insert() to populate.
+        explicit CellSet(CellSetID id);
 
-        void identify(CellSetID id);
+        /// Construct non-empty cell set.
+        CellSet(CellSetID id, const std::vector<int>& cells);
 
         const CellSetID& id() const;
 
         void insert(const int cell);
 
+        using const_iterator = IndexSet::const_iterator;
         const_iterator begin() const;
         const_iterator end()   const;
 

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/CellSetValues.cpp

@@ -1,60 +0,0 @@
-/*
-  Copyright 2016 SINTEF ICT, Applied Mathematics.
-  Copyright 2016 Statoil ASA.
-
-  This file is part of the Open Porous Media Project (OPM).
-
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#if HAVE_CONFIG_H
-#include <config.h>
-#endif // HAVE_CONFIG_H
-
-#include <opm/flowdiagnostics/CellSetValues.hpp>
-
-#include <cassert>
-#include <utility>
-
-namespace Opm {
-namespace FlowDiagnostics {
-
-CellSetValues::CellSetValues(const SizeType initialCapacity)
-{
-    assoc_.reserve(initialCapacity);
-}
-
-void
-CellSetValues::addCellValue(const int    cellIndex,
-                                 const double cellValue)
-{
-    assoc_.push_back(Association{cellIndex, cellValue});
-}
-
-CellSetValues::SizeType
-CellSetValues::cellValueCount() const
-{
-    return assoc_.size();
-}
-
-CellSetValues::CellValue
-CellSetValues::cellValue(const SizeType cellValueIndex) const
-{
-    const auto& a = assoc_[cellValueIndex];
-
-    return { a.index, a.value };
-}
-
-} // namespace FlowDiagnostics
-} // namespace Opm

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/CellSetValues.hpp

@@ -21,54 +21,12 @@
 #ifndef OPM_CELLSETVALUES_HEADER_INCLUDED
 #define OPM_CELLSETVALUES_HEADER_INCLUDED
 
-#include <utility>
+#include <map>
-#include <vector>
 
 namespace Opm {
 namespace FlowDiagnostics {
 
-    class CellSetValues
+    using CellSetValues = std::map<int, double>;
-    {
-    public:
-        using SizeType  = std::vector<int>::size_type;
-        using CellValue = std::pair<int, double>;
-
-        /// Constructor.
-        ///
-        /// @param[in] initialCapacity Number of elements that can be stored
-        ///            in set without reallocation.
-        explicit CellSetValues(const SizeType initialCapacity = 0);
-
-        /// Associate value with specific cell, represented by its index.
-        ///
-        /// @param[in] cellIndex Index of specific cell.
-        ///
-        /// @param[in] cellValue Value associated with cell @p cellIndex.
-        void addCellValue(const int    cellIndex,
-                          const double cellValue);
-
-        /// Retrieve number of elements stored in set.
-        SizeType cellValueCount() const;
-
-        /// Retrieve value association for single set element.
-        ///
-        /// @param[in] cellValueIndex Linear ID of single cell->value
-        ///            association.  Must be in the range @code [0,
-        ///            cellValueCount()) @endcode.
-        ///
-        /// @returns Single association between cell index and numerical
-        /// value.
-        CellValue cellValue(const SizeType cellValueIndex) const;
-
-    private:
-        struct Association
-        {
-            int    index;
-            double value;
-        };
-
-        std::vector<Association> assoc_;
-    };
 
 } // namespace FlowDiagnostics
 } // namespace Opm

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.cpp

@@ -0,0 +1,254 @@
+/*
+  Copyright 2015, 2016 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#include <opm/flowdiagnostics/DerivedQuantities.hpp>
+#include <algorithm>
+#include <numeric>
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+    namespace {
+        /// Helper function for flowCapacityStorageCapacityCurve().
+        template <class InputValues, class ExtractElement>
+        std::vector<double>
+        cumulativeNormalized(const InputValues& input,
+                             ExtractElement&& extract)
+        {
+            // Extract quantity.
+            auto q = std::vector<double>{};
+            q.reserve(input.size() + 1);
+            q.push_back(0.0);
+            for (const auto& e : input) {
+                q.push_back(extract(e));
+            }
+
+            // Accumulate and normalize.
+            std::partial_sum(q.begin(), q.end(), q.begin());
+            const auto t = q.back();
+            for (auto& qi : q) {
+                qi /= t;
+            }
+
+            return q;
+        }
+
+    } // anonymous namespace
+
+
+
+    /// The F-Phi curve.
+    ///
+    /// 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).
+    Graph flowCapacityStorageCapacityCurve(const Toolbox::Forward& injector_solution,
+                                           const Toolbox::Reverse& producer_solution,
+                                           const std::vector<double>& pv)
+    {
+        const auto& ftof = injector_solution.fd.timeOfFlight();
+        const auto& rtof = producer_solution.fd.timeOfFlight();
+        if (pv.size() != ftof.size() || pv.size() != rtof.size()) {
+            throw std::runtime_error("flowCapacityStorageCapacityCurve(): "
+                                     "Input solutions 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());
+
+        auto Phi = cumulativeNormalized(time_and_pv, [](const D2& i) { return i.first; });
+        auto F = cumulativeNormalized(time_and_pv, [](const D2& i) { return i.second / i.first; });
+
+        return Graph{std::move(Phi), std::move(F)};
+    }
+
+
+
+
+    /// The Lorenz coefficient from 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.
+    double lorenzCoefficient(const Graph& flowcap_storagecap_curve)
+    {
+        const auto& storagecap = flowcap_storagecap_curve.first;
+        const auto& flowcap = flowcap_storagecap_curve.second;
+        if (flowcap.size() != storagecap.size()) {
+            throw std::runtime_error("lorenzCoefficient(): Inconsistent sizes in input graph.");
+        }
+        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);
+    }
+
+
+
+
+    /// Compute sweep efficiency versus dimensionless time (pore
+    /// volumes injected).
+    ///
+    /// The sweep efficiency is analogue to 1D displacement using
+    /// the F-Phi curve as flux function.
+    Graph sweepEfficiency(const Graph& flowcap_storagecap_curve)
+    {
+        const auto& storagecap = flowcap_storagecap_curve.first;
+        const auto& flowcap = flowcap_storagecap_curve.second;
+        if (flowcap.size() != storagecap.size()) {
+            throw std::runtime_error("sweepEfficiency(): Inconsistent sizes in input graph.");
+        }
+
+        // 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(tD, Ev);
+    }
+
+
+
+
+
+    /// Compute pore volume associated with an injector-producer pair.
+    double injectorProducerPairVolume(const Toolbox::Forward& injector_solution,
+                                      const Toolbox::Reverse& producer_solution,
+                                      const std::vector<double>& pore_volume,
+                                      const CellSetID& injector,
+                                      const CellSetID& producer)
+    {
+        const auto& inj_tracer = injector_solution.fd.concentration(injector);
+        const auto& prod_tracer = producer_solution.fd.concentration(producer);
+
+        double volume = 0.0;
+        for (const auto& inj_data : inj_tracer) {
+            const int cell = inj_data.first;
+            const auto prod_data = prod_tracer.find(cell);
+            if (prod_data != prod_tracer.end()) {
+                volume += pore_volume[cell] * inj_data.second * prod_data->second;
+            }
+        }
+        return volume;
+    }
+
+
+
+
+
+    namespace {
+
+        // Helper for injectorProducerPairFlux().
+        double pairFlux(const CellSetValues& tracer,
+                        const CellSet& well_cells,
+                        const CellSetValues& inflow_flux,
+                        const bool require_inflow)
+        {
+            double flux = 0.0;
+            for (const int cell : well_cells) {
+                const auto tracer_iter = tracer.find(cell);
+                if (tracer_iter != tracer.end()) {
+                    // Tracer present in cell.
+                    const auto source_iter = inflow_flux.find(cell);
+                    if (source_iter != inflow_flux.end()) {
+                        // Cell has source term.
+                        const double source = source_iter->second;
+                        if ((source > 0.0) == require_inflow) {
+                            // Source term has correct sign.
+                            flux += source * tracer_iter->second;
+                        }
+                    }
+                }
+            }
+            return flux;
+        }
+
+    } // anonymous namespace
+
+
+
+
+    /// Compute fluxes associated with an injector-producer pair.
+    ///
+    /// The first flux returned is the injection flux associated with the given producers,
+    /// (equal to the accumulated product of producer tracer values at the injector cells
+    /// with the injection fluxes), the second is the production flux associated with the
+    /// given injectors. In general, they will only be the same (up to sign) for
+    /// incompressible cases.
+    ///
+    /// Note: since we consider injecting fluxes positive and producing fluxes negative
+    /// (for the inflow_flux), the first returned element will be positive and the second
+    /// will be negative.
+    std::pair<double, double>
+    injectorProducerPairFlux(const Toolbox::Forward& injector_solution,
+                             const Toolbox::Reverse& producer_solution,
+                             const CellSet& injector_cells,
+                             const CellSet& producer_cells,
+                             const CellSetValues& inflow_flux)
+    {
+        const auto& inj_tracer = injector_solution.fd.concentration(injector_cells.id());
+        const auto& prod_tracer = producer_solution.fd.concentration(producer_cells.id());
+        const double inj_flux = pairFlux(prod_tracer, injector_cells, inflow_flux, true);
+        const double prod_flux = pairFlux(inj_tracer, producer_cells, inflow_flux, false);
+        return { inj_flux, prod_flux };
+    }
+
+
+} // namespace FlowDiagnostics
+} // namespace Opm

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.hpp

@@ -0,0 +1,106 @@
+/*
+  Copyright 2016 Statoil ASA.
+  Copyright 2015, 2016 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_DERIVEDQUANTITIES_HEADER_INCLUDED
+#define OPM_DERIVEDQUANTITIES_HEADER_INCLUDED
+
+#include <opm/flowdiagnostics/CellSet.hpp>
+#include <opm/flowdiagnostics/Solution.hpp>
+#include <opm/flowdiagnostics/Toolbox.hpp>
+#include <utility>
+#include <vector>
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+    /// Class used to return graph objects. For a graph g, the
+    /// abscissa (x) values should go in g.first and the ordinate (y)
+    /// values in g.second.
+    using Graph = std::pair< std::vector<double>, std::vector<double> >;
+
+    /// The F-Phi curve.
+    ///
+    /// 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).
+    ///
+    /// Returns F (flow capacity) as a function of Phi (storage capacity),
+    /// that is for the returned Graph g, g.first is Phi and g.second is F.
+    Graph flowCapacityStorageCapacityCurve(const Toolbox::Forward& injector_solution,
+                                           const Toolbox::Reverse& producer_solution,
+                                           const std::vector<double>& pore_volume);
+
+    /// The Lorenz coefficient from 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.
+    double lorenzCoefficient(const Graph& flowcap_storagecap_curve);
+
+    /// Compute sweep efficiency versus dimensionless time (pore
+    /// volumes injected).
+    ///
+    /// The sweep efficiency is analogue to 1D displacement using
+    /// the F-Phi curve as flux function.
+    Graph sweepEfficiency(const Graph& flowcap_storagecap_curve);
+
+    /// Compute pore volume associated with an injector-producer pair.
+    double injectorProducerPairVolume(const Toolbox::Forward& injector_solution,
+                                      const Toolbox::Reverse& producer_solution,
+                                      const std::vector<double>& pore_volume,
+                                      const CellSetID& injector,
+                                      const CellSetID& producer);
+
+    /// Compute fluxes associated with an injector-producer pair.
+    ///
+    /// The first flux returned is the injection flux associated with the given producers,
+    /// (equal to the accumulated product of producer tracer values at the injector cells
+    /// with the injection fluxes), the second is the production flux associated with the
+    /// given injectors. In general, they will only be the same (up to sign) for
+    /// incompressible cases.
+    ///
+    /// Note: since we consider injecting fluxes positive and producing fluxes negative
+    /// (for the inflow_flux), the first returned element will be positive and the second
+    /// will be negative.
+    std::pair<double, double>
+    injectorProducerPairFlux(const Toolbox::Forward& injector_solution,
+                             const Toolbox::Reverse& producer_solution,
+                             const CellSet& injector_cells,
+                             const CellSet& producer_cells,
+                             const CellSetValues& inflow_flux);
+
+
+} // namespace FlowDiagnostics
+} // namespace Opm
+
+
+
+#endif // OPM_DERIVEDQUANTITIES_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/Solution.cpp

@@ -53,17 +53,8 @@ public:
     CellSetValues concentration(const CellSetID& tracer) const;
 
 private:
-    struct CompareCellSetIDs
-    {
-        bool operator()(const CellSetID& x,
-                        const CellSetID& y) const
-        {
-            return x.to_string() < y.to_string();
-        }
-    };
-
     using SolutionMap =
-        std::map<CellSetID, CellSetValues, CompareCellSetIDs>;
+        std::map<CellSetID, CellSetValues>;
 
     GlobalToF   tof_;
     SolutionMap tracerToF_;

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/Toolbox.cpp

@@ -69,8 +69,8 @@ private:
     CellSetValues       only_inflow_flux_;
     CellSetValues       only_outflow_flux_;
 
-    AssembledConnections inj_conn_;
+    AssembledConnections downstream_conn_;
-    AssembledConnections prod_conn_;
+    AssembledConnections upstream_conn_;
     bool conn_built_ = false;
 
     void buildAssembledConnections();
@@ -111,26 +111,13 @@ Toolbox::Impl::assignConnectionFlux(const ConnectionValues& flux)
 void
 Toolbox::Impl::assignInflowFlux(const CellSetValues& inflow_flux)
 {
-    // Count the inflow (>0) fluxes.
+    only_inflow_flux_.clear();
-    const int num_items = inflow_flux.cellValueCount();
+    only_outflow_flux_.clear();
-    int num_inflows = 0;
+    for (const auto& data : inflow_flux) {
-    for (int item = 0; item < num_items; ++item) {
-        if (inflow_flux.cellValue(item).second > 0.0) {
-            ++num_inflows;
-        }
-    }
-
-    // Reserve memory.
-    only_inflow_flux_ = CellSetValues(num_inflows);
-    only_outflow_flux_ = CellSetValues(num_items - num_inflows);
-
-    // Build in- and out-flow structures.
-    for (int item = 0; item < num_items; ++item) {
-        auto data = inflow_flux.cellValue(item);
         if (data.second > 0.0) {
-            only_inflow_flux_.addCellValue(data.first, data.second);
+            only_inflow_flux_[data.first] = data.second;
         } else {
-            only_outflow_flux_.addCellValue(data.first, -data.second);
+            only_outflow_flux_[data.first] = -data.second;
         }
     }
 }
@@ -151,7 +138,7 @@ Toolbox::Impl::injDiag(const std::vector<CellSet>& start_sets)
     using ToF = Solution::TimeOfFlight;
     using Conc = Solution::TracerConcentration;
 
-    TracerTofSolver solver(inj_conn_, pvol_, only_inflow_flux_);
+    TracerTofSolver solver(downstream_conn_, upstream_conn_, pvol_, only_inflow_flux_);
     sol.assignGlobalToF(solver.solveGlobal(start_sets));
 
     for (const auto& start : start_sets) {
@@ -179,7 +166,7 @@ Toolbox::Impl::prodDiag(const std::vector<CellSet>& start_sets)
     using ToF = Solution::TimeOfFlight;
     using Conc = Solution::TracerConcentration;
 
-    TracerTofSolver solver(prod_conn_, pvol_, only_outflow_flux_);
+    TracerTofSolver solver(upstream_conn_, downstream_conn_, pvol_, only_outflow_flux_);
     sol.assignGlobalToF(solver.solveGlobal(start_sets));
 
     for (const auto& start : start_sets) {
@@ -197,8 +184,8 @@ Toolbox::Impl::buildAssembledConnections()
     // Create the data structures needed by the tracer/tof solver.
     const size_t num_connections = g_.numConnections();
     const size_t num_phases = flux_.numPhases();
-    inj_conn_ = AssembledConnections();
+    downstream_conn_ = AssembledConnections();
-    prod_conn_ = AssembledConnections();
+    upstream_conn_ = AssembledConnections();
     for (size_t conn_idx = 0; conn_idx < num_connections; ++conn_idx) {
         auto cells = g_.connection(conn_idx);
         using ConnID = ConnectionValues::ConnID;
@@ -209,16 +196,16 @@ Toolbox::Impl::buildAssembledConnections()
             connection_flux += flux_(ConnID{conn_idx}, PhaseID{phase});
         }
         if (connection_flux > 0.0) {
-            inj_conn_.addConnection(cells.first, cells.second, connection_flux);
+            downstream_conn_.addConnection(cells.first, cells.second, connection_flux);
-            prod_conn_.addConnection(cells.second, cells.first, connection_flux);
+            upstream_conn_.addConnection(cells.second, cells.first, connection_flux);
         } else {
-            inj_conn_.addConnection(cells.second, cells.first, -connection_flux);
+            downstream_conn_.addConnection(cells.second, cells.first, -connection_flux);
-            prod_conn_.addConnection(cells.first, cells.second, -connection_flux);
+            upstream_conn_.addConnection(cells.first, cells.second, -connection_flux);
         }
     }
     const int num_cells = g_.numCells();
-    inj_conn_.compress(num_cells);
+    downstream_conn_.compress(num_cells);
-    prod_conn_.compress(num_cells);
+    upstream_conn_.compress(num_cells);
 
     // Mark as built (until flux changed).
     conn_built_ = true;

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/TracerTofSolver.cpp

@@ -41,9 +41,7 @@ namespace FlowDiagnostics
         std::vector<double> expandSparse(const int n, const CellSetValues& v)
         {
             std::vector<double> r(n, 0.0);
-            const int num_items = v.cellValueCount();
+            for (const auto& data : v) {
-            for (int item = 0; item < num_items; ++item) {
-                auto data = v.cellValue(item);
                 r[data.first] = data.second;
             }
             return r;
@@ -78,9 +76,10 @@ namespace FlowDiagnostics
 
 
     TracerTofSolver::TracerTofSolver(const AssembledConnections& graph,
+                                     const AssembledConnections& reverse_graph,
                                      const std::vector<double>& pore_volumes,
                                      const CellSetValues& source_inflow)
-        : TracerTofSolver(graph, pore_volumes, source_inflow, InOutFluxComputer(graph))
+        : TracerTofSolver(graph, reverse_graph, pore_volumes, source_inflow, InOutFluxComputer(graph))
     {
     }
 
@@ -91,10 +90,12 @@ namespace FlowDiagnostics
     // The InOutFluxComputer is used so that influx_ and outflux_ can be
     // const members of the class.
     TracerTofSolver::TracerTofSolver(const AssembledConnections& graph,
+                                     const AssembledConnections& reverse_graph,
                                      const std::vector<double>& pore_volumes,
                                      const CellSetValues& source_inflow,
                                      InOutFluxComputer&& inout)
         : g_(graph)
+        , g_reverse_(reverse_graph)
         , pv_(pore_volumes)
         , influx_(std::move(inout.influx))
         , outflux_(std::move(inout.outflux))
@@ -138,12 +139,14 @@ namespace FlowDiagnostics
 
         // Return computed time-of-flight.
         CellSetValues local_tof;
+        CellSetValues local_tracer;
         const int num_elements = component_starts_.back();
         for (int element = 0; element < num_elements; ++element) {
             const int cell = sequence_[element];
-            local_tof.addCellValue(cell, tof_[cell]);
+            local_tof[cell] = tof_[cell];
+            local_tracer[cell] = tracer_[cell];
         }
-        return LocalSolution{ local_tof, CellSetValues{} }; // TODO also return tracer
+        return LocalSolution{ std::move(local_tof), std::move(local_tracer) };
     }
 
 
@@ -156,10 +159,10 @@ namespace FlowDiagnostics
         const int num_cells = pv_.size();
         is_start_.clear();
         is_start_.resize(num_cells, 0);
-        upwind_contrib_.clear();
-        upwind_contrib_.resize(num_cells, 0.0);
         tof_.clear();
-        tof_.resize(num_cells, -1e100);
+        tof_.resize(num_cells, 0.0);
+        tracer_.clear();
+        tracer_.resize(num_cells, 0.0);
         num_multicell_ = 0;
         max_size_multicell_ = 0;
         max_iter_multicell_ = 0;
@@ -263,6 +266,11 @@ namespace FlowDiagnostics
                 solveMultiCell(comp_size, &sequence_[component_starts_[comp]]);
             }
         }
+
+        // Threshold time-of-flight values.
+        for (double& t : tof_) {
+            t = std::min(t, max_tof_);
+        }
     }
 
 
@@ -272,28 +280,38 @@ namespace FlowDiagnostics
     void TracerTofSolver::solveSingleCell(const int cell)
     {
         // Compute influx (divisor of tof expression).
-        double source = 2.0 * source_term_[cell];  // Initial tof for well cell equal to half fill time.
+        double source = source_term_[cell];  // Initial tof for well cell equal to fill time.
         if (source == 0.0 && is_start_[cell]) {
             source = std::numeric_limits<double>::infinity(); // Gives 0 tof in start cell.
         }
         const double total_influx = influx_[cell] + source;
 
-        // Compute effective pv (dividend of tof expression).
+        // Cap time-of-flight if time to fill cell is greater than
-        const double eff_pv = pv_[cell] + upwind_contrib_[cell];
+        // max_tof_. Note that cells may still have larger than
-
+        // max_tof_ after solveSingleCell() when including upwind
-        // Compute (capped) tof.
+        // contributions, and those in turn can affect cells
-        if (total_influx < eff_pv / max_tof_) {
+        // downstream (so capping in this method will not produce the
+        // same result). All tofs will finally be capped in solve() as
+        // a post-process. The reason for the somewhat convoluted
+        // behaviour is to match existing MRST results.
+        if (total_influx < pv_[cell] / max_tof_) {
             tof_[cell] = max_tof_;
-        } else {
+            return;
-            tof_[cell] = eff_pv / total_influx;
         }
 
-        // Set contribution for my downwind cells (if any).
+        // Compute upwind contribution.
-        for (const auto& conn : g_.cellNeighbourhood(cell)) {
+        double upwind_tof_contrib = 0.0;
-            const int downwind_cell = conn.neighbour;
+        double upwind_tracer_contrib = 0.0;
+        for (const auto& conn : g_reverse_.cellNeighbourhood(cell)) {
+            const int upwind_cell = conn.neighbour;
             const double flux = conn.weight;
-            upwind_contrib_[downwind_cell] += tof_[cell] * flux;
+            upwind_tof_contrib += tof_[upwind_cell] * flux;
+            upwind_tracer_contrib += tracer_[upwind_cell] * flux;
         }
+
+        // Compute time-of-flight and tracer.
+        tof_[cell] = (pv_[cell] + upwind_tof_contrib) / total_influx;
+        tracer_[cell] = is_start_[cell] ? 1.0 : upwind_tracer_contrib / total_influx;
     }
 
 
@@ -302,9 +320,9 @@ namespace FlowDiagnostics
 
     void TracerTofSolver::solveMultiCell(const int num_cells, const int* cells)
     {
+        // Record some statistics.
         ++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;
@@ -318,7 +336,6 @@ namespace FlowDiagnostics
                 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);
     }

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/TracerTofSolver.hpp

@@ -49,8 +49,10 @@ namespace FlowDiagnostics
     public:
         /// Initialize solver with a given flow graph (a weighted,
         /// directed asyclic graph) containing the out-fluxes from
-        /// each cell, pore volumes and inflow sources (positive).
+        /// each cell, the reverse graph (with in-fluxes from each
+        /// cell), pore volumes and inflow sources (positive).
         TracerTofSolver(const AssembledConnections& graph,
+                        const AssembledConnections& reverse_graph,
                         const std::vector<double>& pore_volumes,
                         const CellSetValues& source_inflow);
 
@@ -77,6 +79,7 @@ namespace FlowDiagnostics
         // --------------  Private data members --------------
 
         const AssembledConnections& g_;
+        const AssembledConnections& g_reverse_;
         const std::vector<double>& pv_;
         const std::vector<double> influx_;
         const std::vector<double> outflux_;
@@ -84,8 +87,8 @@ namespace FlowDiagnostics
         std::vector<char> is_start_; // char to avoid the nasty vector<bool> specialization
         std::vector<int> sequence_;
         std::vector<int> component_starts_;
-        std::vector<double> upwind_contrib_;
         std::vector<double> tof_;
+        std::vector<double> tracer_;
         int num_multicell_ = 0;
         int max_size_multicell_ = 0;
         int max_iter_multicell_ = 0;
@@ -99,6 +102,7 @@ namespace FlowDiagnostics
         // --------------  Private methods --------------
 
         TracerTofSolver(const AssembledConnections& graph,
+                        const AssembledConnections& reverse_graph,
                         const std::vector<double>& pore_volumes,
                         const CellSetValues& source_inflow,
                         InOutFluxComputer&& inout);

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_cellset.cpp

@@ -56,6 +56,12 @@ BOOST_AUTO_TEST_CASE (Construct)
 
         BOOST_CHECK_EQUAL(i.to_string(), name);
     }
+    {
+        const auto i1 = CellSetID("I-1");
+        const auto i2 = CellSetID("I-2");
+        BOOST_CHECK_EQUAL(i1 < i2, true);
+        BOOST_CHECK_EQUAL(i1 < i2, i1.to_string() < i2.to_string());
+    }
 }
 
 BOOST_AUTO_TEST_SUITE_END()
@@ -65,19 +71,10 @@ BOOST_AUTO_TEST_SUITE(CellSetTest)
 
 BOOST_AUTO_TEST_CASE (Constructor)
 {
-    {
-        auto s = CellSet{};
-
-        BOOST_CHECK_EQUAL(s.id().to_string(), "");
-    }
-
     {
         const auto name = std::string("Test-Ctor");
 
-        auto s = CellSet{};
+        auto s = CellSet{CellSetID(name)};
-        {
-            s.identify(CellSetID(name));
-        }
 
         BOOST_CHECK_EQUAL(s.id().to_string(), name);
     }
@@ -85,27 +82,43 @@ BOOST_AUTO_TEST_CASE (Constructor)
 
 BOOST_AUTO_TEST_CASE (AssignCells)
 {
-    auto s = CellSet{};
-
     const auto cells = std::vector<int>
         { 0, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000 };
 
-    for (const auto& cell : cells) {
-        s.insert(cell);
-    }
-
-    auto out = std::vector<int>(s.begin(), s.end());
     {
-        std::sort(out.begin(), out.end());
+        // Using insert() to populate.
+        auto s = CellSet{CellSetID("TestSet")};
+        for (const auto& cell : cells) {
+            s.insert(cell);
+        }
+
+        auto out = std::vector<int>(s.begin(), s.end());
+        {
+            std::sort(out.begin(), out.end());
+        }
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(out  .begin(), out  .end(),
+                                      cells.begin(), cells.end());
     }
 
-    BOOST_CHECK_EQUAL_COLLECTIONS(out  .begin(), out  .end(),
+    {
-                                  cells.begin(), cells.end());
+        // Using direct constructor to populate.
+        auto s = CellSet{CellSetID("TestSet"), cells};
+
+        auto out = std::vector<int>(s.begin(), s.end());
+        {
+            std::sort(out.begin(), out.end());
+        }
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(out  .begin(), out  .end(),
+                                      cells.begin(), cells.end());
+    }
 }
 
+
 BOOST_AUTO_TEST_CASE (Duplicates)
 {
-    auto s = CellSet{};
+    auto s = CellSet{CellSetID("TestSet")};
 
     const auto cells = std::vector<int>
         { 0, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000 };
@@ -125,4 +138,25 @@ BOOST_AUTO_TEST_CASE (Duplicates)
                                   cells.begin(), cells.end());
 }
 
+
+BOOST_AUTO_TEST_CASE (DuplicatesDirectConstruction)
+{
+    const auto cells = std::vector<int>
+        { 0, 100, 100, 100, 2, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000 };
+
+    const auto expected = std::vector<int>
+        { 0, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000 };
+
+    auto s = CellSet{CellSetID("TestSet"), cells};
+
+    auto out = std::vector<int>(s.begin(), s.end());
+    {
+        std::sort(out.begin(), out.end());
+    }
+
+    BOOST_CHECK_EQUAL_COLLECTIONS(out     .begin(), out     .end(),
+                                  expected.begin(), expected.end());
+}
+
+
 BOOST_AUTO_TEST_SUITE_END()

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_cellsetvalues.cpp

@@ -1,89 +0,0 @@
-/*
-  Copyright 2016 SINTEF ICT, Applied Mathematics.
-  Copyright 2016 Statoil ASA.
-
-  This file is part of the Open Porous Media Project (OPM).
-
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#if HAVE_CONFIG_H
-#include <config.h>
-#endif // HAVE_CONFIG_H
-
-#if HAVE_DYNAMIC_BOOST_TEST
-#define BOOST_TEST_DYN_LINK
-#endif
-
-#define NVERBOSE
-
-#define BOOST_TEST_MODULE TEST_CELLSETVALUES
-
-#include <opm/common/utility/platform_dependent/disable_warnings.h>
-#include <boost/test/unit_test.hpp>
-#include <opm/common/utility/platform_dependent/reenable_warnings.h>
-
-#include <opm/flowdiagnostics/CellSetValues.hpp>
-
-using Opm::FlowDiagnostics::CellSetValues;
-
-BOOST_AUTO_TEST_SUITE(CellSet_Values)
-
-BOOST_AUTO_TEST_CASE (Constructor)
-{
-    {
-        CellSetValues s{};
-    }
-
-    {
-        auto s = CellSetValues{ 100 };
-    }
-}
-
-BOOST_AUTO_TEST_CASE (AssignValues)
-{
-    auto s = CellSetValues{ 100 };
-
-    for (decltype(s.cellValueCount())
-             i = 0, n = 100;
-         i < n; ++i)
-    {
-        s.addCellValue(100 - i, i * 10.0);
-    }
-
-    BOOST_CHECK_EQUAL(s.cellValueCount(), 100);
-
-    {
-        const auto a = s.cellValue(0);
-
-        BOOST_CHECK_EQUAL(a.first , 100);
-        BOOST_CHECK_CLOSE(a.second, 0.0, 1.0e-10);
-    }
-
-    {
-        const auto a = s.cellValue(s.cellValueCount() - 1);
-
-        BOOST_CHECK_EQUAL(a.first , 1);
-        BOOST_CHECK_CLOSE(a.second, 990.0, 1.0e-10);
-    }
-
-    {
-        const auto a = s.cellValue(50);
-
-        BOOST_CHECK_EQUAL(a.first , 50);
-        BOOST_CHECK_CLOSE(a.second, 500.0, 1.0e-10);
-    }
-}
-
-BOOST_AUTO_TEST_SUITE_END()

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_derivedquantities.cpp

@@ -0,0 +1,299 @@
+/*
+  Copyright 2016 Statoil ASA.
+
+  This file is part of the Open Porous Media Project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif // HAVE_DYNAMIC_BOOST_TEST
+
+#define NVERBOSE
+
+#define BOOST_TEST_MODULE TEST_DERIVEDQUANTITIES
+
+#include <boost/test/unit_test.hpp>
+
+#include <opm/flowdiagnostics/DerivedQuantities.hpp>
+#include <opm/flowdiagnostics/Toolbox.hpp>
+
+#include <algorithm>
+
+using namespace Opm::FlowDiagnostics;
+
+namespace
+{
+    std::size_t
+    numIntConn(const std::size_t nx,
+               const std::size_t ny)
+    {
+        return (nx - 1)*ny + nx*(ny - 1);
+    }
+
+    std::vector<int>
+    internalConnections(const std::size_t nx,
+                        const std::size_t ny)
+    {
+        auto cellID = [](const std::size_t start,
+                         const std::size_t off)
+        {
+            return static_cast<int>(start + off);
+        };
+
+        auto neighbours = std::vector<int>{};
+        neighbours.reserve(2 * numIntConn(nx, ny));
+
+        // I connections
+        {
+            for (auto j = 0*ny; j < ny; ++j) {
+                const auto start = j * nx;
+
+                for (auto i = 0*nx + 1; i < nx; ++i) {
+                    neighbours.push_back(cellID(start, i - 1));
+                    neighbours.push_back(cellID(start, i - 0));
+                }
+            }
+        }
+
+        // J connections
+        {
+            for (auto j = 0*ny + 1; j < ny; ++j) {
+                const auto start = (j - 1)*nx;
+
+                for (auto i = 0*nx; i < nx; ++i) {
+                    neighbours.push_back(cellID(start, i + 0 ));
+                    neighbours.push_back(cellID(start, i + nx));
+                }
+            }
+        }
+
+        return neighbours;
+    }
+
+    std::vector<double>
+    flowField(const std::vector<double>::size_type n)
+    {
+        return std::vector<double>(n, 0.3);
+    }
+
+} // Namespace anonymous
+
+class Setup
+{
+public:
+    Setup(const std::size_t nx,
+          const std::size_t ny);
+
+    const ConnectivityGraph& connectivity() const;
+    const std::vector<double>&    poreVolume()   const;
+    const ConnectionValues&  flux()         const;
+
+private:
+    ConnectivityGraph g_;
+    std::vector<double>    pvol_;
+    ConnectionValues  flux_;
+};
+
+Setup::Setup(const std::size_t nx,
+             const std::size_t ny)
+    : g_   (nx * ny, internalConnections(nx, ny))
+    , pvol_(g_.numCells(), 0.3)
+    , flux_(ConnectionValues::NumConnections{ g_.numConnections() },
+            ConnectionValues::NumPhases     { 1 })
+{
+    const auto flux = flowField(g_.numConnections());
+
+    using ConnID = ConnectionValues::ConnID;
+
+    const auto phaseID =
+        ConnectionValues::PhaseID{ 0 };
+
+    for (decltype(flux_.numConnections())
+             conn = 0, nconn = flux_.numConnections();
+         conn < nconn; ++conn)
+    {
+        flux_(ConnID{conn}, phaseID) = flux[conn];
+    }
+}
+
+const ConnectivityGraph&
+Setup::connectivity() const
+{
+    return g_;
+}
+
+const std::vector<double>&
+Setup::poreVolume() const
+{
+    return pvol_;
+}
+
+const ConnectionValues&
+Setup::flux() const
+{
+    return flux_;
+}
+
+BOOST_AUTO_TEST_SUITE(Test_DerivedQuantities)
+
+BOOST_AUTO_TEST_CASE (Constructor)
+{
+    const auto cas = Setup(2, 2);
+
+    Toolbox diagTool(cas.connectivity());
+
+    diagTool.assignPoreVolume(cas.poreVolume());
+    diagTool.assignConnectionFlux(cas.flux());
+}
+
+
+namespace {
+
+    template <typename T>
+    void element_is_close(const T& t1, const T& t2)
+    {
+        BOOST_CHECK_CLOSE(t1, t2, 1.0e-10);
+    }
+
+    // using DP = std::pair<double, double>;
+
+    // template<>
+    // void element_is_close<DP>(const DP& p1, const DP& p2)
+    // {
+    //     BOOST_CHECK_CLOSE(p1.first, p2.first, 1.0e-10);
+    //     BOOST_CHECK_CLOSE(p1.second, p2.second, 1.0e-10);
+    // }
+
+    template <class Collection1, class Collection2>
+    void check_is_close(const Collection1& c1, const Collection2& c2)
+    {
+        BOOST_REQUIRE_EQUAL(c1.size(), c2.size());
+
+        if (! c1.empty()) {
+            auto i1 = c1.begin(), e1 = c1.end();
+            auto i2 = c2.begin();
+
+            for (; i1 != e1; ++i1, ++i2) {
+                element_is_close(*i1, *i2);
+            }
+        }
+    }
+
+    template <>
+    void check_is_close<>(const Graph& c1, const Graph& c2)
+    {
+        BOOST_TEST_MESSAGE("Comparing first collections");
+        check_is_close(c1.first, c2.first);
+        BOOST_TEST_MESSAGE("Comparing second collections");
+        check_is_close(c1.second, c2.second);
+    }
+
+} // Namespace Anonymous
+
+
+
+
+
+
+BOOST_AUTO_TEST_CASE (OneDimCase)
+{
+    using namespace Opm::FlowDiagnostics;
+
+    const auto cas = Setup(5, 1);
+    const auto& graph = cas.connectivity();
+    const auto& pv = cas.poreVolume();
+    const auto& flux = cas.flux();
+
+    // Create well in/out flows.
+    CellSetValues wellflow = { {0, 0.3}, {4, -0.3} };
+
+    Toolbox diagTool(graph);
+    diagTool.assignPoreVolume(pv);
+    diagTool.assignConnectionFlux(flux);
+    diagTool.assignInflowFlux(wellflow);
+
+    auto inje = std::vector<CellSet>{CellSet(CellSetID("I-1"), {0})};
+
+    auto prod = std::vector<CellSet>{CellSet(CellSetID("P-1"), {int(graph.numCells()) - 1})};
+
+    {
+        const auto fwd = diagTool.computeInjectionDiagnostics(inje);
+        const auto rev = diagTool.computeProductionDiagnostics(prod);
+
+        BOOST_TEST_MESSAGE("==== F-Phi graph");
+        const Graph expectedFPhi{
+            { 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 },
+            { 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 }
+        };
+        BOOST_CHECK_THROW(flowCapacityStorageCapacityCurve({}, rev, pv), std::runtime_error);
+        BOOST_CHECK_THROW(flowCapacityStorageCapacityCurve(fwd, {}, pv), std::runtime_error);
+        BOOST_CHECK_THROW(flowCapacityStorageCapacityCurve(fwd, rev, {}), std::runtime_error);
+        const auto fcapscap = flowCapacityStorageCapacityCurve(fwd, rev, pv);
+        check_is_close(fcapscap, expectedFPhi);
+
+        BOOST_TEST_MESSAGE("==== Lorenz coefficient");
+        const double expectedLorenz = 0.0;
+        BOOST_CHECK_CLOSE(lorenzCoefficient(fcapscap), expectedLorenz, 1e-10);
+        const Graph wrongGraph {
+            { 0.0, 0.5, 1.0 },
+            { 1.0, 1.0 }
+        };
+        BOOST_CHECK_THROW(lorenzCoefficient(wrongGraph), std::runtime_error);
+        const Graph maxLorenzGraph {
+            { 0.0, 1.0 },
+            { 1.0, 1.0 }
+        };
+        BOOST_CHECK_CLOSE(lorenzCoefficient(maxLorenzGraph), 1.0, 1e-10);
+        const Graph inbetweenLorenzGraph {
+            { 0.0, 0.45, 1.0 },
+            { 0.0, 0.75, 1.0 }
+        };
+        BOOST_CHECK_CLOSE(lorenzCoefficient(inbetweenLorenzGraph), 0.3, 1e-10);
+
+        BOOST_TEST_MESSAGE("==== Sweep efficiency");
+        const Graph expectedSweep{
+            { 0.0, 1.0, 1.0, 1.0, 1.0, 1.0 },
+            { 0.0, 1.0, 1.0, 1.0, 1.0, 1.0 },
+        };
+        BOOST_CHECK_THROW(sweepEfficiency(wrongGraph), std::runtime_error);
+        check_is_close(sweepEfficiency(fcapscap), expectedSweep);
+        const Graph expSweepMax {
+            { 0.0 },
+            { 0.0 }
+        };
+        check_is_close(sweepEfficiency(maxLorenzGraph), expSweepMax);
+        const Graph expSweepInbetween { // Verified against MRST version
+            { 0.0, 0.6, 2.2 },
+            { 0.0, 0.6, 1.0 }
+        };
+        check_is_close(sweepEfficiency(inbetweenLorenzGraph), expSweepInbetween);
+
+        const double expectedVol12 = 1.5;
+        const double vol12 = injectorProducerPairVolume(fwd, rev, pv, CellSetID("I-1"), CellSetID("P-1"));
+        BOOST_CHECK_CLOSE(vol12, expectedVol12, 1e-10);
+
+        const auto pairflux = injectorProducerPairFlux(fwd, rev, inje[0], prod[0], wellflow);
+        BOOST_CHECK_CLOSE(pairflux.first, 0.3, 1e-10);
+        BOOST_CHECK_CLOSE(pairflux.second, -0.3, 1e-10);
+    }
+
+}
+
+BOOST_AUTO_TEST_SUITE_END()

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_flowdiagnosticstool.cpp

@@ -179,24 +179,8 @@ BOOST_AUTO_TEST_CASE (InjectionDiagnostics)
     diagTool.assignPoreVolume(cas.poreVolume());
     diagTool.assignConnectionFlux(cas.flux());
 
-    auto start = std::vector<CellSet>{};
+    auto start = std::vector<CellSet>{ CellSet(CellSetID("I-1"), {0}),
-    {
+                                       CellSet(CellSetID("I-2"), {int(cas.connectivity().numCells()) - 1}) };
-        start.emplace_back();
-
-        auto& s = start.back();
-
-        s.identify(CellSetID("I-1"));
-        s.insert(0);
-    }
-
-    {
-        start.emplace_back();
-
-        auto& s = start.back();
-
-        s.identify(CellSetID("I-2"));
-        s.insert(cas.connectivity().numCells() - 1);
-    }
 
     const auto fwd = diagTool
         .computeInjectionDiagnostics(start);
@@ -232,15 +216,8 @@ BOOST_AUTO_TEST_CASE (InjectionDiagnostics)
         const auto tof = fwd.fd
             .timeOfFlight(CellSetID("I-1"));
 
-        for (decltype(tof.cellValueCount())
+        for (const auto& v : tof) {
-                 i = 0, n = tof.cellValueCount();
+            BOOST_TEST_MESSAGE("ToF[" << v.first << "] = " << v.second);
-             i < n; ++i)
-        {
-            const auto v = tof.cellValue(i);
-
-            BOOST_TEST_MESSAGE("[" << i << "] -> ToF["
-                               << v.first << "] = "
-                               << v.second);
         }
     }
 
@@ -249,25 +226,18 @@ BOOST_AUTO_TEST_CASE (InjectionDiagnostics)
         const auto conc = fwd.fd
             .concentration(CellSetID("I-2"));
 
-        BOOST_TEST_MESSAGE("conc.cellValueCount() = " <<
+        BOOST_TEST_MESSAGE("conc.size() = " <<
-                           conc.cellValueCount());
+                           conc.size());
 
-        for (decltype(conc.cellValueCount())
+        for (const auto& v : conc) {
-                 i = 0, n = conc.cellValueCount();
+            BOOST_TEST_MESSAGE("Conc[" << v.first << "] = " << v.second);
-             i < n; ++i)
-        {
-            const auto v = conc.cellValue(i);
-
-            BOOST_TEST_MESSAGE("[" << i << "] -> Conc["
-                               << v.first << "] = "
-                               << v.second);
         }
+
     }
 }
 
 
 
-
 namespace {
 
     template <class Collection1, class Collection2>
@@ -308,33 +278,17 @@ BOOST_AUTO_TEST_CASE (OneDimCase)
     }
 
     // Create well in/out flows.
-    CellSetValues wellflow;
+    CellSetValues wellflow = { {0, 0.3}, {4, -0.3} };
-    wellflow.addCellValue(0, 0.3);
-    wellflow.addCellValue(4, -0.3);
 
     Toolbox diagTool(graph);
     diagTool.assignPoreVolume(cas.poreVolume());
     diagTool.assignConnectionFlux(flux);
     diagTool.assignInflowFlux(wellflow);
 
-    auto start = std::vector<CellSet>{};
+    const int first_cell = 0;
-    {
+    const int last_cell = cas.connectivity().numCells() - 1;
-        start.emplace_back();
+    auto start = std::vector<CellSet>{ CellSet(CellSetID("I-1"), {first_cell}),
-
+                                       CellSet(CellSetID("I-2"), {last_cell}) };
-        auto& s = start.back();
-
-        s.identify(CellSetID("I-1"));
-        s.insert(0);
-    }
-
-    {
-        start.emplace_back();
-
-        auto& s = start.back();
-
-        s.identify(CellSetID("I-2"));
-        s.insert(cas.connectivity().numCells() - 1);
-    }
 
     const auto fwd = diagTool.computeInjectionDiagnostics(start);
     const auto rev = diagTool.computeProductionDiagnostics(start);
@@ -344,7 +298,7 @@ BOOST_AUTO_TEST_CASE (OneDimCase)
         const auto tof = fwd.fd.timeOfFlight();
 
         BOOST_REQUIRE_EQUAL(tof.size(), cas.connectivity().numCells());
-        std::vector<double> expected = { 0.5, 1.5, 2.5, 3.5, 0.0 };
+        std::vector<double> expected = { 1.0, 2.0, 3.0, 4.0, 0.0 };
         check_is_close(tof, expected);
     }
 
@@ -353,7 +307,7 @@ BOOST_AUTO_TEST_CASE (OneDimCase)
         const auto tof = rev.fd.timeOfFlight();
 
         BOOST_REQUIRE_EQUAL(tof.size(), cas.connectivity().numCells());
-        std::vector<double> expected = { 0.0, 3.5, 2.5, 1.5, 0.5 };
+        std::vector<double> expected = { 0.0, 4.0, 3.0, 2.0, 1.0 };
         check_is_close(tof, expected);
     }
 
@@ -381,15 +335,16 @@ BOOST_AUTO_TEST_CASE (OneDimCase)
         const auto tof = fwd.fd
             .timeOfFlight(CellSetID("I-2"));
 
-        for (decltype(tof.cellValueCount())
-                 i = 0, n = tof.cellValueCount();
-             i < n; ++i)
-        {
-            const auto v = tof.cellValue(i);
 
-            BOOST_TEST_MESSAGE("[" << i << "] -> ToF["
+        std::vector<std::pair<int, double>> expected = { {last_cell, 0.0} };
-                               << v.first << "] = "
+        BOOST_REQUIRE_EQUAL(tof.size(), expected.size());
-                               << v.second);
+
+        int i = 0;
+        for (const auto& v : tof) {
+            BOOST_TEST_MESSAGE("ToF[" << v.first << "] = " << v.second);
+            BOOST_CHECK_EQUAL(v.first, expected[i].first);
+            BOOST_CHECK_CLOSE(v.second, expected[i].second, 1.0e-10);
+            ++i;
         }
     }
 
@@ -398,20 +353,90 @@ BOOST_AUTO_TEST_CASE (OneDimCase)
         const auto conc = fwd.fd
             .concentration(CellSetID("I-2"));
 
-        BOOST_TEST_MESSAGE("conc.cellValueCount() = " <<
+        std::vector<std::pair<int, double>> expected = { {last_cell, 1.0} };
-                           conc.cellValueCount());
+        BOOST_REQUIRE_EQUAL(conc.size(), expected.size());
 
-        for (decltype(conc.cellValueCount())
+        int i = 0;
-                 i = 0, n = conc.cellValueCount();
+        for (const auto& v : conc) {
-             i < n; ++i)
+            BOOST_TEST_MESSAGE("Conc[" << v.first << "] = " << v.second);
-        {
+            BOOST_CHECK_EQUAL(v.first, expected[i].first);
-            const auto v = conc.cellValue(i);
+            BOOST_CHECK_CLOSE(v.second, expected[i].second, 1.0e-10);
+            ++i;
+        }
 
-            BOOST_TEST_MESSAGE("[" << i << "] -> Conc["
+    }
-                               << v.first << "] = "
+
-                               << v.second);
+
+    // Tracer-ToF
+    {
+        const auto tof = fwd.fd
+            .timeOfFlight(CellSetID("I-1"));
+
+        std::vector<double> expected = { 1.0, 2.0, 3.0, 4.0, 5.0 };
+        BOOST_REQUIRE_EQUAL(tof.size(), expected.size());
+
+        for (const auto& v : tof) {
+            BOOST_TEST_MESSAGE("ToF[" << v.first << "] = " << v.second);
+            BOOST_CHECK_CLOSE(v.second, expected[v.first], 1.0e-10);
         }
     }
+
+    // Tracer Concentration
+    {
+        const auto conc = fwd.fd
+            .concentration(CellSetID("I-1"));
+
+        std::vector<double> expected = { 1.0, 1.0, 1.0, 1.0, 1.0 };
+        BOOST_REQUIRE_EQUAL(conc.size(), expected.size());
+
+        for (const auto& v : conc) {
+            BOOST_TEST_MESSAGE("Conc[" << v.first << "] = " << v.second);
+            BOOST_CHECK_CLOSE(v.second, expected[v.first], 1.0e-10);
+        }
+    }
+
+
+    // Add a start point in the middle.
+    const int middle_cell = 2;
+    start.emplace_back(CellSet(CellSetID("Middle"), {middle_cell}));
+
+    const auto fwd2 = diagTool.computeInjectionDiagnostics(start);
+    const auto rev2 = diagTool.computeProductionDiagnostics(start);
+
+    // Tracer-ToF
+    {
+        const auto tof = fwd2.fd
+            .timeOfFlight(CellSetID("Middle"));
+
+        std::vector<std::pair<double, double>> expected = { {2, 0.0}, {3, 1.0}, {4, 2.0} };
+        BOOST_REQUIRE_EQUAL(tof.size(), expected.size());
+
+        int i = 0;
+        for (const auto& v : tof) {
+            BOOST_TEST_MESSAGE("ToF[" << v.first << "] = " << v.second);
+            BOOST_CHECK_EQUAL(v.first, expected[i].first);
+            BOOST_CHECK_CLOSE(v.second, expected[i].second, 1.0e-10);
+            ++i;
+        }
+    }
+
+    // Tracer Concentration
+    {
+        const auto conc = fwd2.fd
+            .concentration(CellSetID("Middle"));
+
+        std::vector<std::pair<double, double>> expected = { {2, 1.0}, {3, 1.0}, {4, 1.0} };
+        BOOST_REQUIRE_EQUAL(conc.size(), expected.size());
+
+        int i = 0;
+        for (const auto& v : conc) {
+            BOOST_TEST_MESSAGE("Conc[" << v.first << "] = " << v.second);
+            BOOST_CHECK_EQUAL(v.first, expected[i].first);
+            BOOST_CHECK_CLOSE(v.second, expected[i].second, 1.0e-10);
+            ++i;
+        }
+    }
+
 }
 
 BOOST_AUTO_TEST_SUITE_END()