#1372 Update flow diagnostics and flow diagnostics applications

ResInsightVersion.cmake

@@ -11,10 +11,10 @@ set(NRLIB_GITHUB_SHA "ba35d4359882f1c6f5e9dc30eb95fe52af50fd6f")
 set(ERT_GITHUB_SHA "06a39878636af0bc52582430ad0431450e51139c") 
 
 # https://github.com/OPM/opm-flowdiagnostics
-set(OPM_FLOWDIAGNOSTICS_SHA "a14dc4ba1302bcc1e0aeb35c5de6b4bd39bce98")
+set(OPM_FLOWDIAGNOSTICS_SHA "2c5fb55db4c4ded49c14161dd16463e1207da049")
 
 # https://github.com/OPM/opm-flowdiagnostics-applications
-set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "1b521494ce9c09a1286dd0a81a1333caa123c680")
+set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "570601718e7197b751bc3cba60c1e5fb7d842842")
 
 # https://github.com/OPM/opm-parser/blob/master/opm/parser/eclipse/Units/Units.hpp
 # This file was moved from opm-core to opm-parser october 2016

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists_files.cmake

@@ -33,6 +33,7 @@ list (APPEND TEST_SOURCE_FILES
         )
 
 list (APPEND EXAMPLE_SOURCE_FILES
+        examples/computeFlowStorageCurve.cpp
         examples/computeLocalSolutions.cpp
         examples/computeToFandTracers.cpp
         examples/computeTracers.cpp

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeFlowStorageCurve.cpp

@@ -0,0 +1,79 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016, 2017 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
+
+#include "exampleSetup.hpp"
+
+#include <opm/flowdiagnostics/DerivedQuantities.hpp>
+#include <numeric>
+
+// Syntax (typical):
+//   computeFlowStorageCurve case=<ecl_case_prefix> step=<report_number>
+int main(int argc, char* argv[])
+try {
+    example::Setup setup(argc, argv);
+    auto& fdTool = setup.toolbox;
+
+    // Solve for forward and reverse time of flight.
+    std::vector<Opm::FlowDiagnostics::CellSet> start;
+    auto fwd = fdTool.computeInjectionDiagnostics(start);
+    auto rev = fdTool.computeProductionDiagnostics(start);
+
+    // Give disconnected cells zero pore volume.
+    std::vector<int> nbcount(setup.graph.numCells(), 0);
+    for (int nb : setup.graph.neighbours()) {
+        if (nb >= 0) {
+            ++nbcount[nb];
+        }
+    }
+    auto pv = setup.graph.poreVolume();
+    for (size_t i = 0; i < pv.size(); ++i) {
+        if (nbcount[i] == 0) {
+            pv[i] = 0.0;
+        }
+    }
+
+    // Cells that have more than 100 times the average pore volume are
+    // probably aquifers, we ignore them (again by setting pore volume
+    // to zero). If this is the correct thing to do or not could
+    // depend on what you want to use the diagnostic for.
+    const double average_pv = std::accumulate(pv.begin(), pv.end(), 0.0) / double(pv.size());
+    for (double& pvval : pv) {
+        if (pvval > 100.0 * average_pv) {
+            pvval = 0.0;
+        }
+    }
+
+    // Compute graph.
+    auto fphi = Opm::FlowDiagnostics::flowCapacityStorageCapacityCurve(fwd, rev, pv);
+
+    // Write it to standard out.
+    std::cout.precision(16);
+    const int sz = fphi.first.size();
+    for (int i = 0; i < sz; ++i) {
+        std::cout << fphi.first[i] << "    " << fphi.second[i] << '\n';
+    }
+}
+catch (const std::exception& e) {
+    std::cerr << "Caught exception: " << e.what() << '\n';
+}

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeLocalSolutions.cpp

@@ -42,9 +42,9 @@ try {
         std::vector<int> completion_cells;
         completion_cells.reserve(well.completions.size());
         for (const auto& completion : well.completions) {
-            const int grid_index = completion.grid_index;
+            const auto& gridName = completion.gridName;
             const auto& ijk = completion.ijk;
-            const int cell_index = setup.graph.activeCell(ijk, grid_index);
+            const int cell_index = setup.graph.activeCell(ijk, gridName);
             if (cell_index >= 0) {
                 completion_cells.push_back(cell_index);
             }

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp

@@ -31,10 +31,10 @@
 
 #include <opm/utility/ECLFluxCalc.hpp>
 #include <opm/utility/ECLGraph.hpp>
+#include <opm/utility/ECLResultData.hpp>
 #include <opm/utility/ECLWellSolution.hpp>
 
 #include <exception>
-#include <iostream>
 #include <sstream>
 #include <stdexcept>
 #include <string>
@@ -80,38 +80,66 @@ namespace example {
         throw std::invalid_argument(os.str());
     }
 
+
+    template <class FluxCalc>
     inline Opm::FlowDiagnostics::ConnectionValues
-    extractFluxField(const Opm::ECLGraph& G, const bool compute_fluxes)
+    extractFluxField(const Opm::ECLGraph& G,
+                     FluxCalc&&           getFlux)
     {
         using ConnVals = Opm::FlowDiagnostics::ConnectionValues;
+
+        const auto actPh = G.activePhases();
+
         auto flux = ConnVals(ConnVals::NumConnections{G.numConnections()},
-                             ConnVals::NumPhases{3});
+                             ConnVals::NumPhases{actPh.size()});
 
         auto phas = ConnVals::PhaseID{0};
 
-        Opm::ECLFluxCalc calc(G);
+        for (const auto& p : actPh) {
+            const auto pflux = getFlux(p);
 
-        const auto phases = { Opm::ECLGraph::PhaseIndex::Aqua   ,
-                              Opm::ECLGraph::PhaseIndex::Liquid ,
-                              Opm::ECLGraph::PhaseIndex::Vapour };
-
-        for (const auto& p : phases)
-        {
-            const auto pflux = compute_fluxes ? calc.flux(p) : G.flux(p);
             if (! pflux.empty()) {
                 assert (pflux.size() == flux.numConnections());
+
                 auto conn = ConnVals::ConnID{0};
                 for (const auto& v : pflux) {
                     flux(conn, phas) = v;
                     conn.id += 1;
                 }
             }
+
             phas.id += 1;
         }
 
         return flux;
     }
 
+    inline Opm::FlowDiagnostics::ConnectionValues
+    extractFluxField(const Opm::ECLGraph&       G,
+                     const Opm::ECLRestartData& rstrt,
+                     const bool                 compute_fluxes)
+    {
+        if (compute_fluxes) {
+            Opm::ECLFluxCalc calc(G);
+
+            auto getFlux = [&calc, &rstrt]
+                (const Opm::ECLGraph::PhaseIndex p)
+            {
+                return calc.flux(rstrt, p);
+            };
+
+            return extractFluxField(G, getFlux);
+        }
+
+        auto getFlux = [&G, &rstrt]
+            (const Opm::ECLGraph::PhaseIndex p)
+        {
+            return G.flux(rstrt, p);
+        };
+
+        return extractFluxField(G, getFlux);
+    }
+
     template <class WellFluxes>
     Opm::FlowDiagnostics::CellSetValues
     extractWellFlows(const Opm::ECLGraph& G,
@@ -120,9 +148,9 @@ namespace example {
         Opm::FlowDiagnostics::CellSetValues inflow;
         for (const auto& well : well_fluxes) {
             for (const auto& completion : well.completions) {
-                const int grid_index = completion.grid_index;
+                const auto& gridName = completion.gridName;
                 const auto& ijk = completion.ijk;
-                const int cell_index = G.activeCell(ijk, grid_index);
+                const int cell_index = G.activeCell(ijk, gridName);
                 if (cell_index >= 0) {
                     // Since inflow is a std::map, if the key was not
                     // already present operator[] will insert a
@@ -142,7 +170,7 @@ namespace example {
 
     struct FilePaths
     {
-        FilePaths(const Opm::parameter::ParameterGroup& param)
+        FilePaths(const Opm::ParameterGroup& param)
         {
             const string casename = param.getDefault<string>("case", "DEFAULT_CASE_NAME");
             grid = param.has("grid") ? param.get<string>("grid")
@@ -164,13 +192,13 @@ namespace example {
 
 
 
-    inline Opm::parameter::ParameterGroup
+    inline Opm::ParameterGroup
     initParam(int argc, char** argv)
     {
         // Obtain parameters from command line (possibly specifying a parameter file).
         const bool verify_commandline_syntax = true;
         const bool parameter_output = false;
-        Opm::parameter::ParameterGroup param(argc, argv, verify_commandline_syntax, parameter_output);
+        Opm::ParameterGroup param(argc, argv, verify_commandline_syntax, parameter_output);
         return param;
     }
 
@@ -180,10 +208,9 @@ namespace example {
     inline Opm::ECLGraph
     initGraph(const FilePaths& file_paths)
     {
-        // Read graph and assign restart file.
+        const auto I = Opm::ECLInitFileData(file_paths.init);
-        auto graph = Opm::ECLGraph::load(file_paths.grid, file_paths.init);
+
-        graph.assignFluxDataSource(file_paths.restart);
+        return Opm::ECLGraph::load(file_paths.grid, I);
-        return graph;
     }
 
 
@@ -209,38 +236,47 @@ namespace example {
     struct Setup
     {
         Setup(int argc, char** argv)
-            : param(initParam(argc, argv))
+            : param          (initParam(argc, argv))
-            , file_paths(param)
+            , file_paths     (param)
-            , graph(initGraph(file_paths))
+            , rstrt          (file_paths.restart)
-            , well_fluxes()
+            , graph          (initGraph(file_paths))
-            , toolbox(initToolbox(graph))
+            , well_fluxes    ()
+            , toolbox        (initToolbox(graph))
             , compute_fluxes_(param.getDefault("compute_fluxes", false))
         {
             const int step = param.getDefault("step", 0);
-            if (!selectReportStep(step)) {
+
+            if (! this->selectReportStep(step)) {
                 std::ostringstream os;
+
                 os << "Report Step " << step
-                   << " is Not Available in Result Set '"
+                   << " is Not Available in Result Set "
-                   << file_paths.grid.stem() << '\'';
+                   << file_paths.grid.stem();
+
                 throw std::domain_error(os.str());
             }
         }
 
         bool selectReportStep(const int step)
         {
-            if (graph.selectReportStep(step)) {
+            if (! rstrt.selectReportStep(step)) {
-                auto wsol = Opm::ECLWellSolution{};
-                well_fluxes = wsol.solution(graph.rawResultData(), graph.numGrids());;
-                toolbox.assignConnectionFlux(extractFluxField(graph, compute_fluxes_));
-                toolbox.assignInflowFlux(extractWellFlows(graph, well_fluxes));
-                return true;
-            } else {
                 return false;
             }
+
+            {
+                auto wsol = Opm::ECLWellSolution{};
+                well_fluxes = wsol.solution(rstrt, graph.activeGrids());
+            }
+
+            toolbox.assignConnectionFlux(extractFluxField(graph, rstrt, compute_fluxes_));
+            toolbox.assignInflowFlux(extractWellFlows(graph, well_fluxes));
+
+            return true;
         }
 
-        Opm::parameter::ParameterGroup param;
+        Opm::ParameterGroup param;
         FilePaths file_paths;
+        Opm::ECLRestartData rstrt;
         Opm::ECLGraph graph;
         std::vector<Opm::ECLWellSolution::WellData> well_fluxes;
         Opm::FlowDiagnostics::Toolbox toolbox;

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/extractFromRestart.cpp

@@ -25,20 +25,21 @@
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/utility/ECLResultData.hpp>
 
+#include <iostream>
 
 // Syntax (typical):
 //   extractFromRestart unrst=<ecl_unrst file> step=<report_number> keyword=<keyword to dump> grid_id=<grid number>
 int main(int argc, char* argv[]) {
     try {
-        Opm::parameter::ParameterGroup param(argc, argv,
+        Opm::ParameterGroup param(argc, argv,
                 /*verify_commandline_syntax=*/ true,
                 /*parameter_output=*/ false);
         const std::string unrst_file = param.get<std::string>("unrst");
         const int report_step = param.getDefault("step", int(0));
-        const int grid_id = param.getDefault("grid_id", int(0));
+        const std::string grid_id = param.getDefault("grid_id", std::string(""));
         const std::string keyword = param.get<std::string>("keyword");
 
-        Opm::ECLResultData restart_file(unrst_file);
+        Opm::ECLRestartData restart_file(unrst_file);
 
         if (!restart_file.selectReportStep(report_step)) {
             std::cerr << "Could not find report step " << report_step << "." << std::endl;

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp

@@ -18,6 +18,8 @@
 */
 
 #include <opm/utility/ECLFluxCalc.hpp>
+#include <opm/utility/ECLResultData.hpp>
+
 #include <opm/parser/eclipse/Units/Units.hpp>
 
 namespace Opm
@@ -34,11 +36,15 @@ namespace Opm
 
 
 
-    std::vector<double> ECLFluxCalc::flux(const PhaseIndex /* phase */) const
+    std::vector<double>
+    ECLFluxCalc::flux(const ECLRestartData& rstrt,
+                      const PhaseIndex /* phase */) const
     {
         // Obtain dynamic data.
         DynamicData dyn_data;
-        dyn_data.pressure = graph_.linearisedCellData("PRESSURE", &ECLUnits::UnitSystem::pressure);
+        dyn_data.pressure = graph_
+            .linearisedCellData(rstrt, "PRESSURE",
+                                &ECLUnits::UnitSystem::pressure);
 
         // Compute fluxes per connection.
         const int num_conn = transmissibility_.size();

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp

@@ -25,6 +25,7 @@
 
 namespace Opm
 {
+    class ECLRestartData;
 
     /// Class for computing connection fluxes in the absence of flux output.
     class ECLFluxCalc
@@ -45,7 +46,9 @@ namespace Opm
         /// \return Flux values corresponding to selected phase.
         ///         Empty if required data is missing.
         ///         Numerical values in SI units (rm^3/s).
-        std::vector<double> flux(const PhaseIndex phase) const;
+        std::vector<double>
+        flux(const ECLRestartData& rstrt,
+             const PhaseIndex      phase) const;
 
     private:
         struct DynamicData

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp

@@ -1,6 +1,5 @@
 /*
-  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016, 2017 Statoil ASA.
-  Copyright 2016 Statoil ASA.
 
   This file is part of the Open Porous Media Project (OPM).
 
@@ -27,6 +26,7 @@
 #include <array>
 #include <cstddef>
 #include <memory>
+#include <string>
 #include <vector>
 
 #include <boost/filesystem.hpp>
@@ -45,8 +45,6 @@ namespace Opm {
     class ECLGraph
     {
     public:
-        using Path = boost::filesystem::path;
-
         /// Enum for indicating requested phase from the flux() method.
         enum class PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
 
@@ -93,17 +91,12 @@ namespace Opm {
         /// \return Fully formed ECLIPSE connection graph with property
         /// associations.
         static ECLGraph
-        load(const Path& grid, const Path& init);
+        load(const boost::filesystem::path& gridFile,
+             const ECLInitFileData&         init);
 
-        /// Assign source object for phase flux calculation.
+        /// Retrieve number of grids in model.
         ///
-        /// \param[in] src Name of ECL restart file, possibly unified, from
+        /// \return The number of LGR grids plus one (the main grid).
-        ///                which next set of phase fluxes should be retrieved.
-        void assignFluxDataSource(const Path& src);
-
-        /// Retrieve number of grids.
-        ///
-        /// \return   The number of LGR grids plus one (the main grid).
         int numGrids() const;
 
         /// Retrieve active cell ID from (I,J,K) tuple in particular grid.
@@ -120,7 +113,7 @@ namespace Opm {
         ///     outside valid range or if the specific cell identified by \p
         ///     ijk and \p gridID is not actually active.
         int activeCell(const std::array<int,3>& ijk,
-                       const int                gridID = 0) const;
+                       const std::string&       gridID = 0) const;
 
         /// Retrieve number of active cells in graph.
         std::size_t numCells() const;
@@ -134,6 +127,11 @@ namespace Opm {
         /// which flux() will return non-zero values if data available.
         const std::vector<PhaseIndex>& activePhases() const;
 
+        /// Retrieve the simulation scenario's set of active grids.
+        ///
+        /// Mostly for canonical lookup of result data in LGRs.
+        const std::vector<std::string>& activeGrids() const;
+
         /// Retrieve neighbourship relations between active cells.
         ///
         /// The \c i-th connection is between active cells \code
@@ -157,24 +155,6 @@ namespace Opm {
         /// \endcode.
         std::vector<double> transmissibility() const;
 
-        /// Restrict dynamic result set data to single report step.
-        ///
-        /// This method must be called before calling either flux() or
-        /// rawResultData().
-        ///
-        /// \param[in] rptstep Selected temporal vector.  Report-step ID.
-        ///
-        /// \return Whether or not dynamic data for the requested report
-        ///    step exists in the underlying result set identified in method
-        ///    assignFluxDataSource().
-        bool selectReportStep(const int rptstep) const;
-
-        /// Access underlying result set.
-        ///
-        /// The result set dynamic data corresponds to the most recent call
-        /// to method selectReportStep().
-        const ::Opm::ECLResultData& rawResultData() const;
-
         /// Retrieve phase flux on all connections defined by \code
         /// neighbours() \endcode.
         ///
@@ -186,7 +166,8 @@ namespace Opm {
         ///    output to the restart file).  Numerical values in SI units
         ///    (rm^3/s).
         std::vector<double>
-        flux(const PhaseIndex phase) const;
+        flux(const ECLRestartData& rstrt,
+             const PhaseIndex      phase) const;
 
         /// Retrieve result set vector from current view (e.g., particular
         /// report step) linearised on active cells.
@@ -196,9 +177,10 @@ namespace Opm {
         /// \param[in] vector Name of result set vector.
         ///
         /// \return Result set vector linearised on active cells.
-        template <typename T>
+        template <typename T, class ResultSet>
         std::vector<T>
-        rawLinearisedCellData(const std::string& vector) const;
+        rawLinearisedCellData(const ResultSet&   rset,
+                              const std::string& vector) const;
 
         /// Convenience type alias for \c UnitSystem PMFs (pointer to member
         /// function).
@@ -211,9 +193,13 @@ namespace Opm {
         /// Typical call:
         /// \code
         ///  const auto press =
-        ///      lCD("PRESSURE", &ECLUnits::UnitSystem::pressure);
+        ///      lCD(rstrt, "PRESSURE", &ECLUnits::UnitSystem::pressure);
         /// \endcode
         ///
+        /// \param[in] rstrt ECL Restart dataset.  It is the responsibility
+        ///    of the caller to ensure that the restart data is correctly
+        ///    positioned on a particular report step.
+        ///
         /// \param[in] vector Name of result set vector.
         ///
         /// \param[in] unit Call-back hook in \c UnitSystem implementation
@@ -223,8 +209,9 @@ namespace Opm {
         /// \return Result set vector linearised on active cells, converted
         ///    to strict SI unit conventions.
         std::vector<double>
-        linearisedCellData(const std::string& vector,
+        linearisedCellData(const ECLRestartData& rstrt,
-                           UnitConvention     unit) const;
+                           const std::string&    vector,
+                           UnitConvention        unit) const;
 
     private:
         /// Implementation class.

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.hpp

@@ -1,6 +1,5 @@
 /*
-  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016, 2017 Statoil ASA.
-  Copyright 2016 Statoil ASA.
 
   This file is part of the Open Porous Media Project (OPM).
 
@@ -18,8 +17,8 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#ifndef OPM_ECLRESTARTDATA_HEADER_INCLUDED
+#ifndef OPM_ECLRESULTDATA_HEADER_INCLUDED
-#define OPM_ECLRESTARTDATA_HEADER_INCLUDED
+#define OPM_ECLRESULTDATA_HEADER_INCLUDED
 
 #include <memory>
 #include <string>
@@ -34,14 +33,17 @@
 /// Forward-declaration of ERT's representation of an ECLIPSE result file.
 ///
 /// This is a hole in the insulation between the interface and the
-/// underlying implementation of class ECLResultData.
+/// underlying implementation of class ECLInitFileData and furthermore
+/// enables constructing wrapper objects from separately parsed result sets.
 extern "C" {
     typedef struct ecl_file_struct ecl_file_type;
 } // extern "C"
 
 namespace Opm {
 
-    /// Representation of an ECLIPSE result-set.
+    class ECLGraph;
+
+    /// Representation of an ECLIPSE Restart result-set.
     ///
     /// This class is aware of the internal structure of ECLIPSE restart
     /// files and may restrict its operation to a single report step.  The
@@ -51,30 +53,37 @@ namespace Opm {
     ///
     /// Note: The client must select a view of the result-set before
     /// accessing any vectors within the set.
-    class ECLResultData
+    class ECLRestartData
     {
     public:
-        using Path = boost::filesystem::path;
-
         /// Default constructor disabled.
-        ECLResultData() = delete;
+        ECLRestartData() = delete;
 
         /// Constructor.
         ///
-        /// \param[in] casePrefix Name or prefix of ECL result data.
+        /// Owning semantics.
-        explicit ECLResultData(const Path& casePrefix);
+        ///
+        /// \param[in] rstrt Name or prefix of ECL result data.
+        explicit ECLRestartData(boost::filesystem::path rstrt);
+
+        /// Constructor
+        ///
+        /// Shared ownership of result set.
+        ///
+        /// \param[in] rstrt ECL restart result set.
+        explicit ECLRestartData(std::shared_ptr<ecl_file_type> rstrt);
 
         /// Copy constructor.
         ///
         /// \param[in] rhs Object from which to construct new instance.
-        ECLResultData(const ECLResultData& rhs);
+        ECLRestartData(const ECLRestartData& rhs);
 
         /// Move constructor.
         ///
         /// \param[in,out] rhs Object from which to construct new instance.
         ///    Its internal implementation will be subsumed into the new
         ///    object.
-        ECLResultData(ECLResultData&& rhs);
+        ECLRestartData(ECLRestartData&& rhs);
 
         /// Assignment operator.
         ///
@@ -82,7 +91,7 @@ namespace Opm {
         ///    instance.
         ///
         /// \return \code *this \endcode.
-        ECLResultData& operator=(const ECLResultData& rhs);
+        ECLRestartData& operator=(const ECLRestartData& rhs);
 
         /// Move assignment operator.
         ///
@@ -91,27 +100,10 @@ namespace Opm {
         ///    instance.
         ///
         /// \return \code *this \endcode.
-        ECLResultData& operator=(ECLResultData&& rhs);
+        ECLRestartData& operator=(ECLRestartData&& rhs);
 
         /// Destructor.
-        ~ECLResultData();
+        ~ECLRestartData();
-
-        /// Access the underlying ERT representation of the result-set.
-        ///
-        /// This is essentially a hole in the interface that is intended to
-        /// support a few very specialised uses.  Handle with care.
-        ///
-        /// \return Handle to underlying ERT representation of result-set.
-        const ecl_file_type* getRawFilePtr() const;
-
-        /// Reset the object's internal view of the result-set to encompass
-        /// the entirety of the underlying file's result vectors.
-        ///
-        /// This is mostly useful in the case of accessing static result
-        /// sets such as INIT files.
-        ///
-        /// \return Whether or not generating the global view succeeded.
-        bool selectGlobalView();
 
         /// Select a result-set view that corresponds to a single report
         /// step.
@@ -123,7 +115,7 @@ namespace Opm {
         /// \return Whether or not selecting the report step succeeded.  The
         ///    typical failure case is the report step not being available
         ///    in the result-set.
-        bool selectReportStep(const int step);
+        bool selectReportStep(const int step) const;
 
         /// Query current result-set view for availability of particular
         /// named result vector in particular enumerated grid.
@@ -132,12 +124,12 @@ namespace Opm {
         ///    availability.
         ///
         /// \param[in] gridID Identity of specific grid for which to query
-        ///    data availability.
+        ///    data availability.  Empty for the main grid.
         ///
         /// \return Whether or not keyword data for the named result vector
         ///    is available in the specific grid.
         bool haveKeywordData(const std::string& vector,
-                             const int          gridID) const;
+                             const std::string& gridID = "") const;
 
         /// Retrieve current result-set view's data values for particular
         /// named result vector in particular enumerated grid.
@@ -159,13 +151,13 @@ namespace Opm {
         ///    keyword data.
         ///
         /// \param[in] gridID Identity of specific grid for which to
-        ///    retrieve keyword data.
+        ///    retrieve keyword data.  Empty for the main grid.
         ///
         /// \return Keyword data values.  Empty if type conversion fails.
         template <typename T>
         std::vector<T>
         keywordData(const std::string& vector,
-                    const int          gridID) const;
+                    const std::string& gridID = "") const;
 
     private:
         class Impl;
@@ -173,6 +165,119 @@ namespace Opm {
         std::unique_ptr<Impl> pImpl_;
     };
 
+    /// Representation of an ECLIPSE Initialization result-set.
+    ///
+    /// This class is aware of the internal structure of ECLIPSE INIT files
+    /// and queries only those objects that pertain to a single grid.
+    class ECLInitFileData
+    {
+    public:
+        ECLInitFileData() = delete;
+
+        /// Constructor.
+        ///
+        /// Construct from filename.  Owning semantics.
+        ///
+        /// \param[in] casePrefix Name or prefix of ECL result data.
+        explicit ECLInitFileData(boost::filesystem::path initFile);
+
+        /// Constructor.
+        ///
+        /// Construct from dataset already input through other means.
+        ///
+        /// Non-owning/shared ownership semantics.
+        explicit ECLInitFileData(std::shared_ptr<ecl_file_type> initFile);
+
+        /// Copy constructor.
+        ///
+        /// \param[in] rhs Object from which to construct new instance.
+        ECLInitFileData(const ECLInitFileData& rhs);
+
+        /// Move constructor.
+        ///
+        /// \param[in,out] rhs Object from which to construct new instance.
+        ///    Its internal implementation will be subsumed into the new
+        ///    object.
+        ECLInitFileData(ECLInitFileData&& rhs);
+
+        /// Assignment operator.
+        ///
+        /// \param[in] rhs Object from which to assign new values to current
+        ///    instance.
+        ///
+        /// \return \code *this \endcode.
+        ECLInitFileData& operator=(const ECLInitFileData& rhs);
+
+        /// Move assignment operator.
+        ///
+        /// \param[in,out] Object from which to assign new instance values.
+        ///    Its internal implementation will be subsumed into this
+        ///    instance.
+        ///
+        /// \return \code *this \endcode.
+        ECLInitFileData& operator=(ECLInitFileData&& rhs);
+
+        /// Destructor.
+        ~ECLInitFileData();
+
+        /// Query current result-set view for availability of particular
+        /// named result vector in particular enumerated grid.
+        ///
+        /// \param[in] vector Named result vector for which to query data
+        ///    availability.
+        ///
+        /// \param[in] gridID Identity of specific grid for which to query
+        ///    data availability.  Empty for the main grid.
+        ///
+        /// \return Whether or not keyword data for the named result vector
+        ///    is available in the specific grid.
+        bool haveKeywordData(const std::string& vector,
+                             const std::string& gridID = "") const;
+
+        /// Retrieve current result-set view's data values for particular
+        /// named result vector in particular enumerated grid.
+        ///
+        /// Will fail (throw an exception of type std::invalid_argument)
+        /// unless the requested keyword data is available in the specific
+        /// grid in the current active view.
+        ///
+        /// \tparam T Element type of return value.  The underlying keyword
+        ///    data will be converted to this type if needed and possible.
+        ///    Note that some type combinations do not make sense.  It is,
+        ///    for instance, not possible to retrieve keyword values of an
+        ///    underlying arithmetic type in the form of a \code
+        ///    std::vector<std::string> \endcode.  Similarly, we cannot
+        ///    access underlying character data through elements of an
+        ///    arithmetic type (e.g., \code std::vector<double> \endcode.)
+        ///
+        /// \param[in] vector Named result vector for which to retrieve
+        ///    keyword data.
+        ///
+        /// \param[in] gridID Identity of specific grid for which to
+        ///    retrieve keyword data.  Empty for the main grid.
+        ///
+        /// \return Keyword data values.  Empty if type conversion fails.
+        template <typename T>
+        std::vector<T>
+        keywordData(const std::string& vector,
+                    const std::string& gridID = "") const;
+
+        // Grant class ECLGraph privileged access to getRawFilePtr().
+        friend class ECLGraph;
+
+    private:
+        class Impl;
+
+        std::unique_ptr<Impl> pImpl_;
+
+        /// Access the underlying ERT representation of the result-set.
+        ///
+        /// This is essentially a hole in the interface that is intended to
+        /// support a few very specialised uses.  Handle with care.
+        ///
+        /// \return Handle to underlying ERT representation of result-set.
+        const ecl_file_type* getRawFilePtr() const;
+    };
 } // namespace Opm
 
-#endif  // OPM_ECLRESTARTDATA_HEADER_INCLUDED
+#endif  // OPM_ECLRESULTDATA_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.cpp

@@ -1,5 +1,4 @@
 /*
-  Copyright 2017 SINTEF ICT, Applied Mathematics.
   Copyright 2017 Statoil ASA.
 
   This file is part of the Open Porous Media Project (OPM).

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.hpp

@@ -1,5 +1,4 @@
 /*
-  Copyright 2017 SINTEF ICT, Applied Mathematics.
   Copyright 2017 Statoil ASA.
 
   This file is part of the Open Porous Media Project (OPM).

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp

@@ -1,6 +1,6 @@
 /*
   Copyright 2016 SINTEF ICT, Applied Mathematics.
-  Copyright 2016 Statoil ASA.
+  Copyright 2016, 2017 Statoil ASA.
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -146,13 +146,16 @@ namespace Opm
 
 
     std::vector<ECLWellSolution::WellData>
-    ECLWellSolution::solution(const ECLResultData& restart,
+    ECLWellSolution::solution(const ECLRestartData& restart,
-                              const int num_grids) const
+                              const std::vector<std::string>& grids) const
     {
+        // Note: this function expects to be called with the correct restart
+        // block--e.g., a report step--selected in the caller.
+
         // Read well data for global grid.
         std::vector<WellData> all_wd{};
-        for (int grid_index = 0; grid_index < num_grids; ++grid_index) {
+        for (const auto& gridName : grids) {
-            std::vector<WellData> wd = readWellData(restart, grid_index);
+            std::vector<WellData> wd = readWellData(restart, gridName);
             // Append to set of all well data.
             all_wd.insert(all_wd.end(), wd.begin(), wd.end());
         }
@@ -163,24 +166,34 @@ namespace Opm
 
 
     std::vector<ECLWellSolution::WellData>
-    ECLWellSolution::readWellData(const ECLResultData& restart, const int grid_index) const
+    ECLWellSolution::readWellData(const ECLRestartData& restart,
+                                  const std::string&    gridName) const
     {
-        // Note: this function is expected to be called in a context
+        // Check if result set provides complete set of well solution data.
-        // where the correct restart block using the ert block mechanisms.
+        if (! (restart.haveKeywordData(ZWEL_KW, gridName) &&
+               restart.haveKeywordData(IWEL_KW, gridName) &&
+               restart.haveKeywordData("XWEL" , gridName) &&
+               restart.haveKeywordData(ICON_KW, gridName) &&
+               restart.haveKeywordData("XCON" , gridName)))
+        {
+            // Not all requisite keywords present in this grid.  Can't
+            // create a well solution.
+            return {};
+        }
 
         // Read header, return if trivial.
-        INTEHEAD ih(restart.keywordData<int>(INTEHEAD_KW, grid_index));
+        INTEHEAD ih(restart.keywordData<int>(INTEHEAD_KW, gridName));
         if (ih.nwell == 0) {
             return {};
         }
         const double qr_unit = resRateUnit(ih.unit);
 
-        // Read necessary keywords.
+        // Load well topology and flow rates.
-        auto zwel = restart.keywordData<std::string>(ZWEL_KW, grid_index);
+        auto zwel = restart.keywordData<std::string>(ZWEL_KW, gridName);
-        auto iwel = restart.keywordData<int>        (IWEL_KW, grid_index);
+        auto iwel = restart.keywordData<int>        (IWEL_KW, gridName);
-        auto xwel = restart.keywordData<double>     ("XWEL" , grid_index);
+        auto xwel = restart.keywordData<double>     ("XWEL" , gridName);
-        auto icon = restart.keywordData<int>        (ICON_KW, grid_index);
+        auto icon = restart.keywordData<int>        (ICON_KW, gridName);
-        auto xcon = restart.keywordData<double>     ("XCON" , grid_index);
+        auto xcon = restart.keywordData<double>     ("XCON" , gridName);
 
         // Create well data.
         std::vector<WellData> wd_vec;
@@ -204,7 +217,7 @@ namespace Opm
                 const int xcon_offset = (well*ih.ncwma + comp_index) * ih.nxcon;
                 WellData::Completion completion;
                 // Note: subtracting 1 from indices (Fortran -> C convention).
-                completion.grid_index = grid_index;
+                completion.gridName = gridName;
                 completion.ijk = { icon[icon_offset + ICON_I_INDEX] - 1,
                                    icon[icon_offset + ICON_J_INDEX] - 1,
                                    icon[icon_offset + ICON_K_INDEX] - 1 };
@@ -212,7 +225,7 @@ namespace Opm
                 completion.reservoir_inflow_rate = -unit::convert::from(xcon[xcon_offset + XCON_QR_INDEX], qr_unit);
                 if (disallow_crossflow_) {
                     // Add completion only if not cross-flowing (injecting producer or producing injector).
-                    if (completion.reservoir_inflow_rate < 0.0 == is_producer) {
+                    if ((completion.reservoir_inflow_rate < 0.0) == is_producer) {
                         wd.completions.push_back(completion);
                     }
                 } else {

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp

@@ -1,6 +1,6 @@
 /*
   Copyright 2016 SINTEF ICT, Applied Mathematics.
-  Copyright 2016 Statoil ASA.
+  Copyright 2016, 2017 Statoil ASA.
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -29,7 +29,7 @@
 namespace Opm
 {
 
-    class ECLResultData;
+    class ECLRestartData;
 
     class ECLWellSolution
     {
@@ -47,7 +47,7 @@ namespace Opm
             bool is_injector_well;
             struct Completion
             {
-                int grid_index;                // 0 for main grid, otherwise LGR grid.
+                std::string gridName;          // empty for main grid, otherwise LGR grid.
                 std::array<int, 3> ijk;        // Cartesian location in grid.
                 double reservoir_inflow_rate;  // Total fluid rate in SI (m^3/s).
             };
@@ -58,8 +58,8 @@ namespace Opm
         ///
         /// Will throw if required data is not available for the
         /// requested step.
-        std::vector<WellData> solution(const ECLResultData& restart,
+        std::vector<WellData> solution(const ECLRestartData& restart,
-                                       const int num_grids) const;
+                                       const std::vector<std::string>& grids) const;
 
     private:
         // Data members.
@@ -67,8 +67,8 @@ namespace Opm
         bool disallow_crossflow_;
 
         // Methods.
-        std::vector<WellData> readWellData(const ECLResultData& restart,
+        std::vector<WellData> readWellData(const ECLRestartData& restart,
-                                           const int grid_index) const;
+                                           const std::string& gridName) const;
     };
 
 

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runAcceptanceTest.cpp

@@ -250,7 +250,7 @@ namespace {
     }
 
     ErrorTolerance
-    testTolerances(const ::Opm::parameter::ParameterGroup& param)
+    testTolerances(const ::Opm::ParameterGroup& param)
     {
         const auto atol = param.getDefault("atol", 1.0e-8);
         const auto rtol = param.getDefault("rtol", 5.0e-12);
@@ -277,7 +277,7 @@ namespace {
     }
 
     ReferenceToF
-    loadReference(const ::Opm::parameter::ParameterGroup& param,
+    loadReference(const ::Opm::ParameterGroup& param,
                   const int                               step,
                   const int                               nDigits)
     {

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runLinearisedCellDataTest.cpp

@@ -377,7 +377,7 @@ namespace {
     }
 
     ErrorTolerance
-    testTolerances(const ::Opm::parameter::ParameterGroup& param)
+    testTolerances(const ::Opm::ParameterGroup& param)
     {
         const auto atol = param.getDefault("atol", 1.0e-8);
         const auto rtol = param.getDefault("rtol", 5.0e-12);
@@ -386,7 +386,7 @@ namespace {
     }
 
     std::vector<std::string>
-    testQuantities(const ::Opm::parameter::ParameterGroup& param)
+    testQuantities(const ::Opm::ParameterGroup& param)
     {
         return StringUtils::VectorValue::
             get<std::string>(param.get<std::string>("quant"));
@@ -411,7 +411,7 @@ namespace {
     }
 
     ReferenceSolution
-    loadReference(const ::Opm::parameter::ParameterGroup& param,
+    loadReference(const ::Opm::ParameterGroup& param,
                   const std::string&                      quant,
                   const int                               step,
                   const int                               nDigits)
@@ -493,7 +493,8 @@ namespace {
 
     std::array<AggregateErrors, 2>
     sampleDifferences(const ::Opm::ECLGraph&                  graph,
-                      const ::Opm::parameter::ParameterGroup& param,
+                      const ::Opm::ECLRestartData&            rstrt,
+                      const ::Opm::ParameterGroup& param,
                       const std::string&                      quant,
                       const std::vector<int>&                 steps)
     {
@@ -510,7 +511,7 @@ namespace {
         auto E = std::array<AggregateErrors, 2>{};
 
         for (const auto& step : steps) {
-            if (! graph.selectReportStep(step)) {
+            if (! rstrt.selectReportStep(step)) {
                 continue;
             }
 
@@ -518,7 +519,7 @@ namespace {
 
             {
                 const auto raw = Calculated {
-                    graph.rawLinearisedCellData<double>(ECLquant)
+                    graph.rawLinearisedCellData<double>(rstrt, ECLquant)
                 };
 
                 computeErrors(Reference{ ref.raw }, raw, E[0]);
@@ -526,7 +527,7 @@ namespace {
 
             {
                 const auto SI = Calculated {
-                    graph.linearisedCellData(ECLquant, unit)
+                    graph.linearisedCellData(rstrt, ECLquant, unit)
                 };
 
                 computeErrors(Reference{ ref.SI }, SI, E[1]);
@@ -561,12 +562,14 @@ try {
     const auto pth = example::FilePaths(prm);
     const auto tol = testTolerances(prm);
 
+    const auto rstrt = ::Opm::ECLRestartData(pth.restart);
     const auto steps = availableReportSteps(pth);
     const auto graph = example::initGraph(pth);
 
     auto all_ok = true;
     for (const auto& quant : testQuantities(prm)) {
-        const auto E = sampleDifferences(graph, prm, quant, steps);
+        const auto E =
+            sampleDifferences(graph, rstrt, prm, quant, steps);
 
         const auto ok =
             everythingFine(E[0], tol) && everythingFine(E[1], tol);

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runTransTest.cpp

@@ -152,7 +152,7 @@ namespace {
     }
 
     ErrorTolerance
-    testTolerances(const ::Opm::parameter::ParameterGroup& param)
+    testTolerances(const ::Opm::ParameterGroup& param)
     {
         const auto atol = param.getDefault("atol", 1.0e-8);
         const auto rtol = param.getDefault("rtol", 5.0e-12);
@@ -161,7 +161,7 @@ namespace {
     }
 
     std::vector<double>
-    loadReference(const ::Opm::parameter::ParameterGroup& param)
+    loadReference(const ::Opm::ParameterGroup& param)
     {
         namespace fs = boost::filesystem;
 
@@ -185,7 +185,7 @@ namespace {
         };
     }
 
-    bool transfieldAcceptable(const ::Opm::parameter::ParameterGroup& param,
+    bool transfieldAcceptable(const ::Opm::ParameterGroup& param,
                               const std::vector<double>&              trans)
     {
         const auto Tref = loadReference(param);

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

@@ -1,5 +1,6 @@
 /*
-  Copyright 2015, 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2015, 2016, 2017 SINTEF ICT, Applied Mathematics.
+  Copyright 2017 Statoil ASA.
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -70,9 +71,26 @@ namespace FlowDiagnostics
                                            const Toolbox::Reverse& producer_solution,
                                            const std::vector<double>& pv)
     {
-        const auto& ftof = injector_solution.fd.timeOfFlight();
+        return flowCapacityStorageCapacityCurve(injector_solution.fd.timeOfFlight(),
-        const auto& rtof = producer_solution.fd.timeOfFlight();
+                                                producer_solution.fd.timeOfFlight(),
-        if (pv.size() != ftof.size() || pv.size() != rtof.size()) {
+                                                pv);
+    }
+
+
+
+
+    /// 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 std::vector<double>& injector_tof,
+                                           const std::vector<double>& producer_tof,
+                                           const std::vector<double>& pv)
+    {
+        if (pv.size() != injector_tof.size() || pv.size() != producer_tof.size()) {
             throw std::runtime_error("flowCapacityStorageCapacityCurve(): "
                                      "Input solutions must have same size.");
         }
@@ -82,12 +100,12 @@ namespace FlowDiagnostics
         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].first = injector_tof[ii] + producer_tof[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 Phi = cumulativeNormalized(time_and_pv, [](const D2& i) { return i.second; });
         auto F = cumulativeNormalized(time_and_pv, [](const D2& i) { return i.second / i.first; });
 
         return Graph{std::move(Phi), std::move(F)};

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

@@ -50,6 +50,14 @@ namespace FlowDiagnostics
                                            const Toolbox::Reverse& producer_solution,
                                            const std::vector<double>& pore_volume);
 
+    /// This overload gets the injector and producer time-of-flight
+    /// directly instead of extracting it from the solution
+    /// objects. It otherwise behaves identically to the other
+    /// overload.
+    Graph flowCapacityStorageCapacityCurve(const std::vector<double>& injector_tof,
+                                           const std::vector<double>& producer_tof,
+                                           const std::vector<double>& pore_volume);
+
     /// The Lorenz coefficient from the F-Phi curve.
     ///
     /// The Lorenz coefficient is a measure of heterogeneity. It

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

@@ -1,5 +1,6 @@
 /*
   Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
 
   This file is part of the Open Porous Media project (OPM).
 

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

@@ -1,6 +1,6 @@
 /*
   Copyright 2016 SINTEF ICT, Applied Mathematics.
-  Copyright 2016 Statoil ASA.
+  Copyright 2016, 2017 Statoil ASA.
 
   This file is part of the Open Porous Media project (OPM).
 

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/utility/graph/AssembledConnectionsIteration.hpp

@@ -1,5 +1,6 @@
 /*
   Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
 
   This file is part of the Open Porous Media project (OPM).
 

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

@@ -247,6 +247,8 @@ BOOST_AUTO_TEST_CASE (OneDimCase)
         BOOST_CHECK_THROW(flowCapacityStorageCapacityCurve(fwd, rev, {}), std::runtime_error);
         const auto fcapscap = flowCapacityStorageCapacityCurve(fwd, rev, pv);
         check_is_close(fcapscap, expectedFPhi);
+        const auto fcapscap2 = flowCapacityStorageCapacityCurve(fwd.fd.timeOfFlight(), rev.fd.timeOfFlight(), pv);
+        check_is_close(fcapscap2, expectedFPhi);
 
         BOOST_TEST_MESSAGE("==== Lorenz coefficient");
         const double expectedLorenz = 0.0;
@@ -296,4 +298,28 @@ BOOST_AUTO_TEST_CASE (OneDimCase)
 
 }
 
+
+
+
+
+BOOST_AUTO_TEST_CASE (GeneralCase)
+{
+    BOOST_TEST_MESSAGE("==== F-Phi graph");
+
+    std::vector<double> pv { 1.0, 2.0, 1.0 };
+    std::vector<double> ftof { 0.0, 2.0, 1.0 };
+    std::vector<double> rtof { 1.0, 2.0, 0.0 };
+    const Graph expectedFPhi{
+        { 0.0, 0.25, 0.5, 1.0 },
+        { 0.0, 0.4, 0.8, 1.0 }
+    };
+    const auto fcapscap = flowCapacityStorageCapacityCurve(ftof, rtof, pv);
+    check_is_close(fcapscap, expectedFPhi);
+
+    BOOST_TEST_MESSAGE("==== Lorenz coefficient");
+    const double expectedLorenz = 0.3;
+    BOOST_CHECK_CLOSE(lorenzCoefficient(fcapscap), expectedLorenz, 1e-10);
+}
+
+
 BOOST_AUTO_TEST_SUITE_END()