From 3c07eafab29ad56af5e9d58f736b51b4510fc839 Mon Sep 17 00:00:00 2001 From: Magne Sjaastad Date: Fri, 12 May 2017 12:17:11 +0200 Subject: [PATCH] #1372 Update flow diagnostics and flow diagnostics applications --- ResInsightVersion.cmake | 4 +- .../CMakeLists_files.cmake | 1 + .../examples/computeFlowStorageCurve.cpp | 79 ++ .../examples/computeLocalSolutions.cpp | 4 +- .../examples/exampleSetup.hpp | 108 +- .../examples/extractFromRestart.cpp | 7 +- .../opm/utility/ECLFluxCalc.cpp | 10 +- .../opm/utility/ECLFluxCalc.hpp | 5 +- .../opm/utility/ECLGraph.cpp | 547 +++++---- .../opm/utility/ECLGraph.hpp | 63 +- .../opm/utility/ECLResultData.cpp | 1031 ++++++++++++++--- .../opm/utility/ECLResultData.hpp | 185 ++- .../opm/utility/ECLUnitHandling.cpp | 1 - .../opm/utility/ECLUnitHandling.hpp | 1 - .../opm/utility/ECLWellSolution.cpp | 47 +- .../opm/utility/ECLWellSolution.hpp | 14 +- .../tests/runAcceptanceTest.cpp | 4 +- .../tests/runLinearisedCellDataTest.cpp | 19 +- .../tests/runTransTest.cpp | 6 +- .../opm/flowdiagnostics/DerivedQuantities.cpp | 30 +- .../opm/flowdiagnostics/DerivedQuantities.hpp | 8 + .../opm/flowdiagnostics/Solution.cpp | 1 + .../opm/flowdiagnostics/TracerTofSolver.cpp | 2 +- .../graph/AssembledConnectionsIteration.hpp | 1 + .../tests/test_derivedquantities.cpp | 26 + 25 files changed, 1674 insertions(+), 530 deletions(-) create mode 100644 ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeFlowStorageCurve.cpp diff --git a/ResInsightVersion.cmake b/ResInsightVersion.cmake index 2f00b8f1e8..7d3d69bfb6 100644 --- a/ResInsightVersion.cmake +++ b/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 diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists_files.cmake b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists_files.cmake index 568ddf04b3..e30cdd67a1 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists_files.cmake +++ b/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 diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeFlowStorageCurve.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeFlowStorageCurve.cpp new file mode 100644 index 0000000000..f29df16e8b --- /dev/null +++ b/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 . +*/ + +#if HAVE_CONFIG_H +#include +#endif + +#include "exampleSetup.hpp" + +#include +#include + +// Syntax (typical): +// computeFlowStorageCurve case= step= +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 start; + auto fwd = fdTool.computeInjectionDiagnostics(start); + auto rev = fdTool.computeProductionDiagnostics(start); + + // Give disconnected cells zero pore volume. + std::vector 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'; +} diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeLocalSolutions.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeLocalSolutions.cpp index 8b5e2b7152..50a2a5598d 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeLocalSolutions.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeLocalSolutions.cpp @@ -42,9 +42,9 @@ try { std::vector 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); } diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp index 5d2f64b4c0..a536c2730e 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp @@ -31,10 +31,10 @@ #include #include +#include #include #include -#include #include #include #include @@ -80,38 +80,66 @@ namespace example { throw std::invalid_argument(os.str()); } + + template 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 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("case", "DEFAULT_CASE_NAME"); grid = param.has("grid") ? param.get("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. - auto graph = Opm::ECLGraph::load(file_paths.grid, file_paths.init); - graph.assignFluxDataSource(file_paths.restart); - return graph; + const auto I = Opm::ECLInitFileData(file_paths.init); + + return Opm::ECLGraph::load(file_paths.grid, I); } @@ -209,38 +236,47 @@ namespace example { struct Setup { Setup(int argc, char** argv) - : param(initParam(argc, argv)) - , file_paths(param) - , graph(initGraph(file_paths)) - , well_fluxes() - , toolbox(initToolbox(graph)) + : param (initParam(argc, argv)) + , file_paths (param) + , rstrt (file_paths.restart) + , graph (initGraph(file_paths)) + , 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 '" - << file_paths.grid.stem() << '\''; + << " is Not Available in Result Set " + << file_paths.grid.stem(); + throw std::domain_error(os.str()); } } bool selectReportStep(const int step) { - if (graph.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 { + if (! rstrt.selectReportStep(step)) { 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 well_fluxes; Opm::FlowDiagnostics::Toolbox toolbox; diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/extractFromRestart.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/extractFromRestart.cpp index 4d38286871..060fba34f7 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/extractFromRestart.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/extractFromRestart.cpp @@ -25,20 +25,21 @@ #include #include +#include // Syntax (typical): // extractFromRestart unrst= step= keyword= grid_id= 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("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("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; diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp index 91a6e3d3f0..1e597a4064 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp @@ -18,6 +18,8 @@ */ #include +#include + #include namespace Opm @@ -34,11 +36,15 @@ namespace Opm - std::vector ECLFluxCalc::flux(const PhaseIndex /* phase */) const + std::vector + 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(); diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp index 340b1974b0..3f230ab23f 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp +++ b/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 flux(const PhaseIndex phase) const; + std::vector + flux(const ECLRestartData& rstrt, + const PhaseIndex phase) const; private: struct DynamicData diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp index d06645c93d..0e3ee28afd 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp @@ -1,6 +1,5 @@ /* - 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). @@ -39,6 +38,7 @@ #include #include #include +#include #include #include @@ -86,6 +86,17 @@ namespace { const ecl_grid_type* getGrid(const ecl_grid_type* G, const int gridID); + /// Retrieve grid name. + /// + /// \param[in] ERT Grid representation. + /// + /// \param[in] gridID Numeric index of requested grid. Zero for the + /// main grid or positive for one of the LGRs. + /// + /// \return Name of grid \p G. Empty for the main grid. + std::string + getGridName(const ecl_grid_type* G, const int gridID); + /// Extract Cartesian dimensions of an ECL grid. /// /// \param[in] G ERT grid instance corresponding to the model's main @@ -99,14 +110,20 @@ namespace { /// Access unit conventions pertaining to single grid in result set. /// + /// \tparam ResultSet Type representing a result set. Must + /// implement methods \code haveKeywordData() \endcode and \code + /// keywordData<>() \endcode. Typically \code Opm::ECLRestartData + /// \endcode or \code Opm::ECLInitFileData \endcode. + /// /// \param[in] rset Result set. /// - /// \param[in] grid_ID Numerical ID of grid. Non-negative. Zero - /// for the main grid and positive for LGRs. + /// \param[in] gridID ID (name) of particular grid. Empty for the + /// main grid. /// - /// \return Unit system convention for \p grid_ID in result set. - auto getUnitSystem(const ::Opm::ECLResultData& rset, - const int grid_ID) + /// \return Unit system convention for \p gridID in result set. + template + auto getUnitSystem(const ResultSet& rset, + const std::string& gridID) -> decltype(::Opm::ECLUnits::createUnitSystem(0)); /// Retrieve global pore-volume vector from INIT source. @@ -117,15 +134,15 @@ namespace { /// /// \param[in] init ERT representation of INIT source. /// - /// \param[in] grid_ID Numerical ID of grid. Non-negative. Zero - /// for the main grid and positive in the case of an LGR. + /// \param[in] gridID ID (name) of grid. Empty for the main grid + /// and non-empty in the case of an LGR. /// /// \return Vector of pore-volumes for all global cells of \p G in /// SI unit conventions (rm^3). std::vector - getPVolVector(const ecl_grid_type* G, - const ::Opm::ECLResultData& init, - const int grid_ID = 0); + getPVolVector(const ecl_grid_type* G, + const ::Opm::ECLInitFileData& init, + const std::string& gridID = ""); /// Extract non-neighbouring connections from ECLIPSE model /// @@ -137,9 +154,10 @@ namespace { /// \return Model's non-neighbouring connections, including those /// between main and local grids. std::vector - loadNNC(const ecl_grid_type* G, - const ::Opm::ECLResultData& init); + loadNNC(const ecl_grid_type* G, + const ecl_file_type* init); + /// Cartesian connections in a model grid. class CartesianGridData { public: @@ -153,15 +171,17 @@ namespace { /// /// \param[in] gridID Numeric identifier of this grid. Zero for /// main grid, positive for LGRs. - CartesianGridData(const ecl_grid_type* G, - const ::Opm::ECLResultData& init, - const int gridID); + CartesianGridData(const ecl_grid_type* G, + const ::Opm::ECLInitFileData& init, + const int gridID); /// Retrieve non-negative numeric ID of grid instance. /// /// \return Constructor's \c gridID parameter. int gridID() const; + const std::string& gridName() const; + /// Retrieve number of active cells in graph. std::size_t numCells() const; @@ -218,14 +238,15 @@ namespace { /// /// \return Numerical values of result set vector, relative to /// global cell numbering of this grid. + template std::vector - cellData(const ::Opm::ECLResultData& src, - const std::string& vector) const; + cellData(const ResultSet& rset, + const std::string& vector) const; - template + template std::vector - activeCellData(const ::Opm::ECLResultData& src, - const std::string& vector) const; + activeCellData(const ResultSet& rset, + const std::string& vector) const; /// Retrieve values of result set vector for all Cartesian /// connections in grid. @@ -238,9 +259,9 @@ namespace { /// \return Numerical values of result set vector attributed to /// all of the grid's Cartesian connections. std::vector - connectionData(const ::Opm::ECLResultData& src, - const std::string& vector, - const double unit) const; + connectionData(const ::Opm::ECLRestartData& rstrt, + const std::string& vector, + const double unit) const; private: /// Facility for deriving Cartesian neighbourship in a grid @@ -424,6 +445,9 @@ namespace { /// than zero for LGRs. const int gridID_; + /// Grid name. Mostly for querying properties in local grids. + const std::string gridName_; + /// Map results from active to global cells. CartesianCells cells_; @@ -445,14 +469,19 @@ namespace { /// Predicate for whether or not a particular result vector is /// defined on the grid's cells. /// - /// \param[in] src Result set. + /// \tparam ResultSet Representation of an ECLIPSE result set. + /// Typically \code Opm::ECLInitFileData \endcode or \code + /// Opm::ECLRestartData \endcode. + /// + /// \param[in] rset Result set. /// /// \param[in] vector Name of result vector. /// - /// \return Whether or not \p vector is defined on model's - /// cells and part of the result set \p src. - bool haveCellData(const ::Opm::ECLResultData& src, - const std::string& keyword) const; + /// \return Whether or not \p vector is defined on model's cells + /// and part of the result set \p src. + template + bool haveCellData(const ResultSet& rset, + const std::string& keyword) const; /// Predicate for whether or not a particular result vector is /// defined on the grid's Cartesian connections. @@ -462,10 +491,10 @@ namespace { /// \param[in] vector Prefix of result vector name. /// /// \return Whether or not all vectors formed by \p vector plus - /// known directional suffixes are defined on model's cells and - /// part of the result set \p src. - bool haveConnData(const ::Opm::ECLResultData& src, - const std::string& keyword) const; + /// known directional suffixes are defined on model's cells + /// and part of the result set \p src. + bool haveConnData(const ::Opm::ECLRestartData& rstrt, + const std::string& keyword) const; /// Append directional cell data to global collection of /// connection data identified by vector name prefix. @@ -480,7 +509,7 @@ namespace { /// input, collection of values corresponding to any previous /// directions (preserved), and on output additionally contains /// the data corresponding to the Cartesian direction \p d. - void connectionData(const ::Opm::ECLResultData& src, + void connectionData(const ::Opm::ECLRestartData& rstrt, const CartesianCells::Direction d, const std::string& vector, const double unit, @@ -510,7 +539,7 @@ namespace { /// /// \param[in] init Internalised void deriveNeighbours(const std::vector& gcells, - const ::Opm::ECLResultData& init, + const ::Opm::ECLInitFileData& init, const CartesianCells::Direction d); }; } // namespace ECL @@ -535,10 +564,20 @@ ECL::getGrid(const ecl_grid_type* G, const int gridID) return ecl_grid_iget_lgr(G, gridID - 1); } +std::string +ECL::getGridName(const ecl_grid_type* G, const int gridID) +{ + if (gridID == ECL_GRID_MAINGRID_LGR_NR) { + return ""; // Empty for main grid. + } + + return ecl_grid_get_name(G); +} + std::vector -ECL::getPVolVector(const ecl_grid_type* G, - const ::Opm::ECLResultData& init, - const int gridID) +ECL::getPVolVector(const ecl_grid_type* G, + const ::Opm::ECLInitFileData& init, + const std::string& gridID) { auto make_szt = [](const int i) { @@ -578,8 +617,8 @@ ECL::loadCase(const boost::filesystem::path& grid) if (! G) { std::ostringstream os; - os << "Failed to load ECL Grid from '" - << grid.generic_string() << '\''; + os << "Failed to load ECL Grid from " + << grid.generic_string(); throw std::invalid_argument(os.str()); } @@ -600,21 +639,40 @@ ECL::cartesianDimensions(const ecl_grid_type* G) make_szt(ecl_grid_get_nz(G)) } }; } -auto ECL::getUnitSystem(const ::Opm::ECLResultData& rset, - const int grid_ID) +template +auto ECL::getUnitSystem(const ResultSet& rset, + const std::string& grid_ID) -> decltype(::Opm::ECLUnits::createUnitSystem(0)) { assert (rset.haveKeywordData(INTEHEAD_KW, grid_ID) && "Result Set Does Not Provide Grid Header"); - const auto ih = rset.keywordData(INTEHEAD_KW, grid_ID); + const auto ih = rset.template keywordData(INTEHEAD_KW, grid_ID); - return ::Opm::ECLUnits::createUnitSystem(ih[INTEHEAD_UNIT_INDEX]); + const auto usys = ih[INTEHEAD_UNIT_INDEX]; + const auto validUsys = (usys >= 1) && (usys <= 4); + + if (! validUsys) { + if (! grid_ID.empty()) { + // Unity system not provided in local grid. Fall back to + // querying the main grid instead. + + const auto mainGrid = std::string{ "" }; + + return getUnitSystem(rset, mainGrid); + } + + throw std::out_of_range { + "No Valid Unit System Key in Result-Set" + }; + } + + return ::Opm::ECLUnits::createUnitSystem(usys); } std::vector -ECL::loadNNC(const ecl_grid_type* G, - const ::Opm::ECLResultData& init) +ECL::loadNNC(const ecl_grid_type* G, + const ecl_file_type* init) { auto make_szt = [](const int n) { @@ -628,7 +686,7 @@ ECL::loadNNC(const ecl_grid_type* G, if (numNNC > 0) { nncData.resize(numNNC); - ecl_nnc_export(G, init.getRawFilePtr(), nncData.data()); + ecl_nnc_export(G, init, nncData.data()); } return nncData; @@ -853,7 +911,7 @@ ECL::CartesianGridData::CartesianCells::numActiveCells() const { return this->rsMap_.subset.size(); } - + std::size_t ECL::CartesianGridData:: CartesianCells::globIdx(const IndexTuple& ijk) const @@ -890,11 +948,12 @@ CartesianCells::ind2sub(const std::size_t globalCell) const // ====================================================================== ECL::CartesianGridData:: -CartesianGridData(const ecl_grid_type* G, - const ::Opm::ECLResultData& init, - const int gridID) - : gridID_(gridID) - , cells_ (G, ::ECL::getPVolVector(G, init, gridID_)) +CartesianGridData(const ecl_grid_type* G, + const ::Opm::ECLInitFileData& init, + const int gridID) + : gridID_ (gridID) + , gridName_(::ECL::getGridName(G, gridID)) + , cells_ (G, ::ECL::getPVolVector(G, init, gridName_)) { { using VT = DirectionSuffix::value_type; @@ -923,6 +982,12 @@ int ECL::CartesianGridData::gridID() const return this->gridID_; } +const std::string& +ECL::CartesianGridData::gridName() const +{ + return this->gridName_; +} + std::size_t ECL::CartesianGridData::numCells() const { @@ -973,48 +1038,51 @@ ECL::CartesianGridData::isSubdivided(const int cellID) const return this->cells_.isSubdivided(cellID); } +template std::vector ECL::CartesianGridData:: -cellData(const ::Opm::ECLResultData& src, - const std::string& vector) const +cellData(const ResultSet& rset, + const std::string& vector) const { - if (! this->haveCellData(src, vector)) { + if (! this->haveCellData(rset, vector)) { return {}; } - auto x = src.keywordData(vector, this->gridID_); + const auto x = + rset.template keywordData(vector, this->gridName()); return this->cells_.scatterToGlobal(x); } namespace { namespace ECL { - template + template std::vector - CartesianGridData::activeCellData(const ::Opm::ECLResultData& src, - const std::string& vector) const + CartesianGridData::activeCellData(const ResultSet& rset, + const std::string& vector) const { - if (! this->haveCellData(src, vector)) { + if (! this->haveCellData(rset, vector)) { return {}; } - auto x = src.keywordData(vector, this->gridID_); + auto x = rset.template keywordData(vector, this->gridName()); return this->cells_.gatherToActive(std::move(x)); } }} // namespace Anonymous::ECL +template bool ECL::CartesianGridData:: -haveCellData(const ::Opm::ECLResultData& src, - const std::string& vector) const +haveCellData(const ResultSet& rset, + const std::string& vector) const { - return src.haveKeywordData(vector, this->gridID_); + return rset.template haveKeywordData(vector, this->gridName()); } bool ECL::CartesianGridData:: -haveConnData(const ::Opm::ECLResultData& src, - const std::string& vector) const +haveConnData(const ::Opm::ECLRestartData& rstrt, + const std::string& vector) const { auto have_data = true; @@ -1024,7 +1092,7 @@ haveConnData(const ::Opm::ECLResultData& src, { const auto vname = this->vectorName(vector, d); - have_data = this->haveCellData(src, vname); + have_data = this->haveCellData(rstrt, vname); if (! have_data) { break; } } @@ -1034,11 +1102,11 @@ haveConnData(const ::Opm::ECLResultData& src, std::vector ECL::CartesianGridData:: -connectionData(const ::Opm::ECLResultData& src, - const std::string& vector, - const double unit) const +connectionData(const ::Opm::ECLRestartData& rstrt, + const std::string& vector, + const double unit) const { - if (! this->haveConnData(src, vector)) { + if (! this->haveConnData(rstrt, vector)) { return {}; } @@ -1048,7 +1116,9 @@ connectionData(const ::Opm::ECLResultData& src, CartesianCells::Direction::J , CartesianCells::Direction::K }) { - this->connectionData(src, d, this->vectorName(vector, d), unit, x); + const auto vname = this->vectorName(vector, d); + + this->connectionData(rstrt, d, vname, unit, x); } return x; @@ -1056,13 +1126,13 @@ connectionData(const ::Opm::ECLResultData& src, void ECL::CartesianGridData:: -connectionData(const ::Opm::ECLResultData& src, +connectionData(const ::Opm::ECLRestartData& rstrt, const CartesianCells::Direction d, const std::string& vector, const double unit, std::vector& x) const { - const auto v = this->cellData(src, vector); + const auto v = this->cellData(rstrt, vector); const auto& cells = this->outCell_.find(d); @@ -1090,7 +1160,7 @@ vectorName(const std::string& vector, void ECL::CartesianGridData:: deriveNeighbours(const std::vector& gcells, - const ::Opm::ECLResultData& init, + const ::Opm::ECLInitFileData& init, const CartesianCells::Direction d) { auto tran = std::string{"TRAN"}; @@ -1112,12 +1182,12 @@ deriveNeighbours(const std::vector& gcells, throw std::invalid_argument("Input direction must be (I,J,K)"); } - const auto& T = init.haveKeywordData(tran, this->gridID_) + const auto& T = init.haveKeywordData(tran, this->gridName()) ? this->cellData(init, tran) : std::vector(this->cells_.numGlobalCells(), 1.0); const auto trans_unit = - ECL::getUnitSystem(init, this->gridID_)->transmissibility(); + ECL::getUnitSystem(init, this->gridName())->transmissibility(); auto SI_trans = [trans_unit](const double trans) { @@ -1166,7 +1236,7 @@ public: /// \param[in] grid Name or prefix of ECL grid (i.e., .GRID or /// .EGRID) file. /// - /// \param[in] init Name of ECL INIT file corresponding to \p grid + /// \param[in] init ECL INIT result set corresponding to \p grid /// input. Assumed to provide at least a complete set /// of pore-volume values (i.e., for all global cells /// defined in the \p grid). @@ -1174,13 +1244,8 @@ public: /// If available in the INIT file, the constructor will /// also leverage the transmissibility data when /// constructing the active cell neighbourship table. - Impl(const Path& grid, const Path& init); - - /// Assign source object for phase flux calculation. - /// - /// \param[in] src Name of ECL restart file, possibly unified, from - /// which next set of phase fluxes should be retrieved. - void assignDataSource(const Path& src); + Impl(const boost::filesystem::path& grid, + const ECLInitFileData& init); /// Retrieve number of grids. /// @@ -1200,7 +1265,7 @@ public: /// ijk from specified grid. Negative one (-1) if (I,J,K) outside /// valid range or if the specific cell identified by \p ijk and \p /// gridID is not actually active. - int activeCell(const int gridID, + int activeCell(const std::string& gridID, const std::array& ijk) const; /// Retrieve number of active cells in graph. @@ -1215,6 +1280,11 @@ public: /// flux() may return non-zero values. const std::vector& activePhases() const; + /// Retrieve the simulation scenario's set of active grids. + /// + /// Mostly for canonical lookup of result data in LGRs. + const std::vector& activeGrids() const; + /// Retrieve neighbourship relations between active cells. /// /// The \c i-th connection is between active cells \code @@ -1238,24 +1308,6 @@ public: /// \endcode. std::vector 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 - /// assignDataSource(). - 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. /// @@ -1269,15 +1321,56 @@ public: /// reported due to report frequencies or no flux values are output at /// all). std::vector - flux(const PhaseIndex phase) const; + flux(const ECLRestartData& rstrt, + const PhaseIndex phase) const; - template + /// Retrieve result set vector from current view linearised on active + /// cells. + /// + /// \tparam T Element type of result set vector. + /// + /// \tparam ResultSet Implementation of an ECL Result Set. Typically + /// \code Opm::ECLRestartData \endcode or \code Opm::ECLInitFileData + /// \endcode. + /// + /// \param[in] rset ECL Result set. Typically an instance of \code + /// Opm::ECLRestartData \endcode or \code Opm::ECLInitFileData + /// \endcode. + /// + /// \param[in] vector Name of result set vector. + /// + /// \return Result set vector linearised on active cells. + template std::vector - rawLinearisedCellData(const std::string& vector) const; + rawLinearisedCellData(const ResultSet& rset, + const std::string& vector) const; + /// Retrieve floating-point result set vector from current view + /// (e.g., particular report step) linearised on active cells and + /// converted to strict SI unit conventions. + /// + /// Typical call: + /// \code + /// const auto press = + /// 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 + /// that enables converting the raw result data to strict SI unit + /// conventions. Hook is called for each grid in the result set. + /// + /// \return Result set vector linearised on active cells, converted + /// to strict SI unit conventions. std::vector - linearisedCellData(const std::string& vector, - UnitConvention unit) const; + linearisedCellData(const ECLRestartData& rstrt, + const std::string& vector, + UnitConvention unit) const; private: /// Collection of non-Cartesian neighbourship relations attributed to a @@ -1518,8 +1611,10 @@ private: /// the simulation run. std::vector activePhases_; - /// Current result set. - std::unique_ptr src_; + /// Set of active grids in result set. + std::vector activeGrids_; + + std::unordered_map gridID_; /// Extract explicit non-neighbouring connections from ECL output. /// @@ -1528,15 +1623,13 @@ private: /// \param[in] G ERT Grid representation. /// /// \param[in] init ERT representation of INIT source. - /// - /// \param[in] coll Backing data for neighbourship extraction. - void defineNNCs(const ecl_grid_type* G, - const ::Opm::ECLResultData& init); + void defineNNCs(const ecl_grid_type* G, + const ECLInitFileData& init); /// Extract scenario's set of active phases. /// /// Writes to activePhases_. - void defineActivePhases(const ::Opm::ECLResultData& init); + void defineActivePhases(const ::Opm::ECLInitFileData& init); /// Compute ECL vector basename for particular phase flux. /// @@ -1554,6 +1647,8 @@ private: /// \tparam[in] GetFluxUnit Type of function object for computing the /// grid-dependent flux unit. /// + /// \param[in] rstrt ECL Restart data result set. + /// /// \param[in] vector Result set vector prefix. Typically computed by /// method flowVector(). /// @@ -1572,9 +1667,10 @@ private: /// contains those values that correspond to the non-neighbouring /// connections (appended onto \p flux). template - void fluxNNC(const std::string& vector, - GetFluxUnit&& fluxUnit, - std::vector& flux) const; + void fluxNNC(const ECLRestartData& rstrt, + const std::string& vector, + GetFluxUnit&& fluxUnit, + std::vector& flux) const; }; // ====================================================================== @@ -1763,12 +1859,10 @@ isViable(const std::vector& grids, // ====================================================================== -Opm::ECLGraph::Impl::Impl(const Path& grid, const Path& init) +Opm::ECLGraph::Impl::Impl(const boost::filesystem::path& grid, + const ECLInitFileData& init) { const auto G = ECL::loadCase(grid); - auto I = ::Opm::ECLResultData{ init }; - - I.selectGlobalView(); const auto numGrids = ECL::numGrids(G.get()); @@ -1779,41 +1873,42 @@ Opm::ECLGraph::Impl::Impl(const Path& grid, const Path& init) for (auto gridID = 0*numGrids; gridID < numGrids; ++gridID) { this->grid_.emplace_back(ECL::getGrid(G.get(), gridID), - I, gridID); + init, gridID); this->activeOffset_.push_back(this->activeOffset_.back() + this->grid_.back().numCells()); + + this->activeGrids_.push_back(this->grid_.back().gridName()); + + this->gridID_[this->activeGrids_.back()] = gridID; } - this->defineNNCs(G.get(), I); - this->defineActivePhases(I); -} - -void -Opm::ECLGraph::Impl::assignDataSource(const Path& src) -{ - this->src_.reset(new ECLResultData(src)); + this->defineNNCs(G.get(), init); + this->defineActivePhases(init); } int -Opm::ECLGraph::Impl:: -numGrids() const +Opm::ECLGraph::Impl::numGrids() const { return grid_.size(); } int Opm::ECLGraph::Impl:: -activeCell(const int gridID, +activeCell(const std::string& gridID, const std::array& ijk) const { - const auto gIdx = - static_castgrid_.size())>(gridID); - - if (gIdx >= this->grid_.size()) { + const auto gID = this->gridID_.find(gridID); + if (gID == std::end(this->gridID_)) { return -1; } + const auto gIdx = + static_castgrid_.size())>(gID->second); + + assert ((gIdx <= this->grid_.size()) && + "Logic Error in ECLGraph::Impl::Impl()"); + const auto& grid = this->grid_[gIdx]; const auto active = grid.activeCell(ijk[0], ijk[1], ijk[2]); @@ -1851,6 +1946,12 @@ Opm::ECLGraph::Impl::activePhases() const return this->activePhases_; } +const std::vector& +Opm::ECLGraph::Impl::activeGrids() const +{ + return this->activeGrids_; +} + std::vector Opm::ECLGraph::Impl::neighbours() const { @@ -1925,24 +2026,13 @@ Opm::ECLGraph::Impl::transmissibility() const return trans; } -const ::Opm::ECLResultData& -Opm::ECLGraph::Impl::rawResultData() const -{ - return *this->src_; -} - -bool Opm::ECLGraph::Impl::selectReportStep(const int rptstep) const -{ - return this->src_->selectReportStep(rptstep); -} - std::vector -Opm::ECLGraph::Impl:: -flux(const PhaseIndex phase) const +Opm::ECLGraph::Impl::flux(const ECLRestartData& rstrt, + const PhaseIndex phase) const { - auto fluxUnit = [this](const int gridID) + auto fluxUnit = [&rstrt](const std::string& gridID) { - return ::ECL::getUnitSystem(*this->src_, gridID)->reservoirRate(); + return ::ECL::getUnitSystem(rstrt, gridID)->reservoirRate(); }; const auto vector = this->flowVector(phase); @@ -1956,7 +2046,7 @@ flux(const PhaseIndex phase) const for (const auto& G : this->grid_) { const auto& q = - G.connectionData(*this->src_, vector, fluxUnit(G.gridID())); + G.connectionData(rstrt, vector, fluxUnit(G.gridName())); if (q.empty()) { // Flux vector invalid unless all grids provide this result @@ -1971,7 +2061,7 @@ flux(const PhaseIndex phase) const // Model includes non-neighbouring connections such as faults and/or // local grid refinement. Extract pertinent flux values for these // connections. - this->fluxNNC(vector, std::move(fluxUnit), v); + this->fluxNNC(rstrt, vector, std::move(fluxUnit), v); } if (v.size() < totconn) { @@ -1985,14 +2075,15 @@ flux(const PhaseIndex phase) const namespace Opm { - template + template std::vector - ECLGraph::Impl::rawLinearisedCellData(const std::string& vector) const + ECLGraph::Impl::rawLinearisedCellData(const ResultSet& rset, + const std::string& vector) const { auto x = std::vector{}; x.reserve(this->numCells()); for (const auto& G : this->grid_) { - const auto xi = G.activeCellData(*this->src_, vector); + const auto xi = G.activeCellData(rset, vector); x.insert(x.end(), std::begin(xi), std::end(xi)); } @@ -2006,19 +2097,20 @@ namespace Opm { } // namespace Opm std::vector -Opm::ECLGraph::Impl::linearisedCellData(const std::string& vector, - UnitConvention unit) const +Opm::ECLGraph::Impl::linearisedCellData(const ECLRestartData& rstrt, + const std::string& vector, + UnitConvention unit) const { auto x = std::vector{}; x.reserve(this->numCells()); for (const auto& G : this->grid_) { - const auto xi = G.activeCellData(*this->src_, vector); + const auto xi = G.activeCellData(rstrt, vector); if (xi.empty()) { continue; } // Note: Compensate for incrementing Grid ID above. const auto usys = - ECL::getUnitSystem(*this->src_, G.gridID()); + ECL::getUnitSystem(rstrt, G.gridName()); // Note: 'usys' (generally, std::unique_ptr<>) does not support // regular PMF syntax (i.e., operator->*()). @@ -2040,25 +2132,28 @@ Opm::ECLGraph::Impl::linearisedCellData(const std::string& vector, } void -Opm::ECLGraph::Impl::defineNNCs(const ecl_grid_type* G, - const ::Opm::ECLResultData& init) +Opm::ECLGraph::Impl::defineNNCs(const ecl_grid_type* G, + const ECLInitFileData& init) { // Assume all transmissibilites in the result set follow the same unit // conventions. - const auto trans_unit = - ECL::getUnitSystem(init, 0)->transmissibility(); + const auto gridID = std::string{ "" }; // Empty in main grid. - for (const auto& nnc : ECL::loadNNC(G, init)) { + const auto trans_unit = + ECL::getUnitSystem(init, gridID)->transmissibility(); + + for (const auto& nnc : ECL::loadNNC(G, init.getRawFilePtr())) { this->nnc_.add(this->grid_, this->activeOffset_, trans_unit, nnc); } } template void -Opm::ECLGraph::Impl::fluxNNC(const std::string& vector, - GetFluxUnit&& fluxUnit, - std::vector& flux) const +Opm::ECLGraph::Impl::fluxNNC(const ECLRestartData& rstrt, + const std::string& vector, + GetFluxUnit&& fluxUnit, + std::vector& flux) const { auto v = std::vector(this->nnc_.numConnections(), 0.0); auto assigned = std::vector(v.size(), false); @@ -2068,28 +2163,31 @@ Opm::ECLGraph::Impl::fluxNNC(const std::string& vector, const auto fluxID = rel.makeKeyword(vector); for (const auto& G : this->grid_) { - const auto gridID = G.gridID(); + const auto& gridName = G.gridName(); const auto& iset = - rel.indexSet().getGridCollection(gridID); + rel.indexSet().getGridCollection(G.gridID()); - if (iset.empty()) { - // No NNCs for this category in this grid. Skip. + if (iset.empty() || + ! rstrt.haveKeywordData(fluxID, gridName)) + { + // No NNCs for this category in this grid or corresponding + // flux vector does not exist. Skip. continue; } - // Note: Method name is confusing, but does actually do what we - // want here. - const auto q = G.cellData(*this->src_, fluxID); + const auto q = rstrt.keywordData(fluxID, gridName); if (q.empty()) { - // No flux data for this category in this grid. Skip. + // Empty flux data for this category in this grid. Not + // really supposed to happen if the above check fires, but + // skip this (category,gridID) pair nonetheless. continue; } - const auto flux_unit = fluxUnit(gridID); + const auto flux_unit = fluxUnit(gridName); - // Data fully available for (category,gridID). Assign + // Data fully available for (category,gridName). Assign // approriate subset of NNC flux vector. for (const auto& ix : iset) { assert (ix.neighIdx < v.size()); @@ -2118,10 +2216,12 @@ Opm::ECLGraph::Impl::fluxNNC(const std::string& vector, void Opm::ECLGraph::Impl:: -defineActivePhases(const ::Opm::ECLResultData& init) +defineActivePhases(const ::Opm::ECLInitFileData& init) { + const auto gridID = std::string{ "" }; // Empty in main grid. + const auto ih = - init.keywordData(INTEHEAD_KW, ECL_GRID_MAINGRID_LGR_NR); + init.keywordData(INTEHEAD_KW, gridID); const auto phaseMask = static_cast(ih[INTEHEAD_PHASE_INDEX]); @@ -2143,8 +2243,7 @@ defineActivePhases(const ::Opm::ECLResultData& init) } std::string -Opm::ECLGraph::Impl:: -flowVector(const PhaseIndex phase) const +Opm::ECLGraph::Impl::flowVector(const PhaseIndex phase) const { const auto vector = std::string("FLR"); // Flow-rate, reservoir @@ -2193,40 +2292,32 @@ Opm::ECLGraph::operator=(ECLGraph&& rhs) } Opm::ECLGraph -Opm::ECLGraph::load(const Path& grid, const Path& init) +Opm::ECLGraph::load(const boost::filesystem::path& grid, + const ECLInitFileData& init) { auto pImpl = ImplPtr{new Impl(grid, init)}; return { std::move(pImpl) }; } -int -Opm::ECLGraph::numGrids() const +int Opm::ECLGraph::numGrids() const { return this->pImpl_->numGrids(); } int Opm::ECLGraph::activeCell(const std::array& ijk, - const int gridID) const + const std::string& gridID) const { return this->pImpl_->activeCell(gridID, ijk); } -void -Opm::ECLGraph::assignFluxDataSource(const Path& src) -{ - this->pImpl_->assignDataSource(src); -} - -std::size_t -Opm::ECLGraph::numCells() const +std::size_t Opm::ECLGraph::numCells() const { return this->pImpl_->numCells(); } -std::size_t -Opm::ECLGraph::numConnections() const +std::size_t Opm::ECLGraph::numConnections() const { return this->pImpl_->numConnections(); } @@ -2237,6 +2328,12 @@ Opm::ECLGraph::activePhases() const return this->pImpl_->activePhases(); } +const std::vector& +Opm::ECLGraph::activeGrids() const +{ + return this->pImpl_->activeGrids(); +} + std::vector Opm::ECLGraph::neighbours() const { return this->pImpl_->neighbours(); @@ -2252,45 +2349,47 @@ std::vector Opm::ECLGraph::transmissibility() const return this->pImpl_->transmissibility(); } -bool Opm::ECLGraph::selectReportStep(const int rptstep) const -{ - return this->pImpl_->selectReportStep(rptstep); -} - -const Opm::ECLResultData& -Opm::ECLGraph::rawResultData() const -{ - return this->pImpl_->rawResultData(); -} - std::vector -Opm::ECLGraph:: -flux(const PhaseIndex phase) const +Opm::ECLGraph::flux(const ECLRestartData& rstrt, + const PhaseIndex phase) const { - return this->pImpl_->flux(phase); + return this->pImpl_->flux(rstrt, phase); } namespace Opm { - template + template std::vector - ECLGraph::rawLinearisedCellData(const std::string& vector) const + ECLGraph::rawLinearisedCellData(const ResultSet& rset, + const std::string& vector) const { - return this->pImpl_->rawLinearisedCellData(vector); + return this->pImpl_->rawLinearisedCellData(rset, vector); } - // Explicit instantiations for the element types we care about. + // Explicit instantiations of method rawLinearisedCellData() for the + // element and result set types we care about. template std::vector - ECLGraph::rawLinearisedCellData(const std::string& vector) const; + ECLGraph::rawLinearisedCellData(const ECLInitFileData& rset, + const std::string& vector) const; + + template std::vector + ECLGraph::rawLinearisedCellData(const ECLRestartData& rset, + const std::string& vector) const; template std::vector - ECLGraph::rawLinearisedCellData(const std::string& vector) const; + ECLGraph::rawLinearisedCellData(const ECLInitFileData& rset, + const std::string& vector) const; + + template std::vector + ECLGraph::rawLinearisedCellData(const ECLRestartData& rset, + const std::string& vector) const; } // namespace Opm std::vector -Opm::ECLGraph::linearisedCellData(const std::string& vector, - UnitConvention unit) const +Opm::ECLGraph::linearisedCellData(const ECLRestartData& rstrt, + const std::string& vector, + UnitConvention unit) const { - return this->pImpl_->linearisedCellData(vector, unit); + return this->pImpl_->linearisedCellData(rstrt, vector, unit); } diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp index 62204e5722..ff60cf67dc 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp @@ -1,6 +1,5 @@ /* - 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). @@ -27,6 +26,7 @@ #include #include #include +#include #include #include @@ -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 - /// 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). + /// \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& 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& activePhases() const; + /// Retrieve the simulation scenario's set of active grids. + /// + /// Mostly for canonical lookup of result data in LGRs. + const std::vector& 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 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 - 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 + template std::vector - 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 - linearisedCellData(const std::string& vector, - UnitConvention unit) const; + linearisedCellData(const ECLRestartData& rstrt, + const std::string& vector, + UnitConvention unit) const; private: /// Implementation class. diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp index 9677d7b180..2add62a4e6 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp @@ -1,6 +1,5 @@ /* - 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). @@ -28,9 +27,15 @@ #include #include #include +#include +#include +#include #include #include +#include +#include #include +#include #include #include @@ -68,9 +73,16 @@ namespace { } namespace ECLImpl { - using FilePtr = ::ERT::ert_unique_ptr; + using FilePtr = std::shared_ptr; namespace Details { + + inline std::string + firstBlockKeyword(const ecl_file_view_type* block) + { + return ecl_kw_get_header(ecl_file_view_iget_kw(block, 0)); + } + /// Convert vector of elements from one element type to another. /// /// \tparam Input Element type of input collection. @@ -378,11 +390,122 @@ namespace { i < nkw; ++i) { result.emplace_back(ecl_kw_iget_char_ptr(kw, i)); + + // Trim trailing white-space. + auto& s = result.back(); + + const auto e = s.find_last_not_of(" \t"); + if (e != std::string::npos) { + s = s.substr(0, e + 1); + } } return result; } }; + + /// Translate grid names to (local) numeric IDs within a + /// section/view of an ECL result set. + /// + /// Provides a caching mechanism to accelerate repeated lookup. + class GridIDCache + { + public: + /// Constructor + /// + /// \param[in] block View relative to which to interpret grid + /// names. + explicit GridIDCache(const ecl_file_view_type* block); + + /// Get integral grid ID of particular grid relative to cache's + /// view. + /// + /// \param[in] gridName Name of particular grid. Empty for the + /// main grid (grid ID 0). + /// + /// \return Numeric grid ID of \p gridName. Negative if \p + /// gridName does not correspond to a grdi in the cache's + /// view. + int getGridID(const std::string& gridName) const; + + private: + /// View into ECL result set. + const ecl_file_view_type* block_; + + /// Cache of name->ID map. + mutable std::unordered_map cache_; + }; + + /// Partition INIT file into sections + class InitFileSections + { + public: + explicit InitFileSections(const ecl_file_type* init); + + struct Section { + Section(const ecl_file_view_type* blk) + : block (blk) + , first_kw(Details::firstBlockKeyword(block)) + {} + + const ecl_file_view_type* block; + std::string first_kw; + }; + + const std::vector
& sections() const + { + return this->sect_; + } + + using SectionID = std::vector
::size_type; + + SectionID numSections() const + { + return this->sections().size(); + } + + const Section& operator[](const SectionID i) const + { + assert ((i < this->numSections()) && "Internal Error"); + + return this->sect_[i]; + } + + private: + const ecl_file_view_type* init_; + + std::vector
sect_; + }; + + template + std::vector getKeywordData(const ecl_kw_type* kw) + { + // Whether or not caller requests a vector. + const auto makeStringVector = + typename std::is_same::type{}; + + switch (getKeywordElementType(kw)) { + case ECL_CHAR_TYPE: + return GetKeywordData:: + as(kw, makeStringVector); + + case ECL_INT_TYPE: + return GetKeywordData:: + as(kw, makeStringVector); + + case ECL_FLOAT_TYPE: + return GetKeywordData:: + as(kw, makeStringVector); + + case ECL_DOUBLE_TYPE: + return GetKeywordData:: + as(kw, makeStringVector); + + default: + // No operator exists for this type. Return empty. + return {}; + } + } } // namespace ECLImpl /// Predicate for whether or not a particular path represents a regular @@ -418,18 +541,22 @@ namespace { /// /// \param[in] file Filename or casename prefix. /// + /// \param[in] ext Set of possible filename extensions. + /// /// \return Filesystem element corresponding to result-set. Either the - /// input \p file itself or, in the case of a casename prefix, the - /// path to a restart file (unified format only). + /// input \p file itself or, in the case of a filename prefix, the + /// first possible match amongst the set generated by the prefix and + /// the input extensions. boost::filesystem::path - deriveResultPath(boost::filesystem::path file) + deriveResultPath(boost::filesystem::path file, + const std::vector& ext) { if (isFile(file)) { return file; } - for (const auto* ext : { ".UNRST", ".FUNRST" }) { - file.replace_extension(ext); + for (const auto& e : ext) { + file.replace_extension(e); if (isFile(file)) { return file; @@ -446,30 +573,72 @@ namespace { throw std::invalid_argument(os.str()); } - /// Load result-set represented by filesystem element. + /// Derive filesystem element from prefix or filename. + /// + /// Pass-through if the input already is a regular file in order to + /// support accessing result-sets other than restart files (e.g., the + /// INIT vectors). /// /// Fails (throws an exception of type \code std::invalid_argument - /// \endcode) if the input filesystem element does not represent a valid - /// result-set resource or if the resource could not be opened (e.g., - /// due to lack of permission). + /// \endcode) if no valid filesystem element can be derived from the + /// input argument /// - /// \param[in] rset Fileystem element representing a result-set. + /// \param[in] file Filename or casename prefix. /// - /// \return Accessor handle of result-set. - ECLImpl::FilePtr openResultSet(const boost::filesystem::path& rset) + /// \return Filesystem element corresponding to result-set. Either the + /// input \p file itself or, in the case of a casename prefix, the + /// path to a restart file (unified format only). + boost::filesystem::path + deriveRestartPath(boost::filesystem::path file) + { + return deriveResultPath(std::move(file), { ".UNRST", ".FUNRST" }); + } + + /// Derive filesystem element from prefix or filename. + /// + /// Pass-through if the input already is a regular file in order to + /// support accessing result-sets other than restart files (e.g., the + /// INIT vectors). + /// + /// Fails (throws an exception of type \code std::invalid_argument + /// \endcode) if no valid filesystem element can be derived from the + /// input argument + /// + /// \param[in] file Filename or casename prefix. + /// + /// \return Filesystem element corresponding to result-set. Either the + /// input \p file itself or, in the case of a casename prefix, the + /// path to an INIT file (unformatted or formatted). + boost::filesystem::path + deriveInitPath(boost::filesystem::path file) + { + return deriveResultPath(std::move(file), { ".INIT", ".FINIT" }); + } + + /// Open ECL result set from pathname. + /// + /// Fails (throws an exception of type \code std::invalid_argument + /// \endcode) if the input argument does not refer to a valid filesystem + /// element. + /// + /// \param[in] file Filename. + /// + /// \return Open stream corresponding to result-set. + ECLImpl::FilePtr openResultSet(const boost::filesystem::path& fname) { // Read-only, keep open between requests const auto open_flags = 0; - auto F = ECLImpl::FilePtr{ - ecl_file_open(rset.generic_string().c_str(), open_flags) + auto F = ECLImpl::FilePtr { + ecl_file_open(fname.generic_string().c_str(), open_flags), + ecl_file_close }; if (! F) { std::ostringstream os; - os << "Failed to load ECL Result object from '" - << rset.generic_string() << '\''; + os << "Failed to load ECL Result object from " + << fname.generic_string(); throw std::invalid_argument(os.str()); } @@ -491,19 +660,160 @@ namespace { auto* globView = ecl_file_get_global_view(const_cast(file)); - return ecl_kw_get_header(ecl_file_view_iget_kw(globView, 0)); + return ECLImpl::Details::firstBlockKeyword(globView); + } + + std::string paddedGridName(const std::string& gridName) + { + if (gridName.empty()) { + return gridName; + } + + std::ostringstream os; + + os << std::setw(8) << std::left << gridName; + + return os.str(); } } // namespace Anonymous -/// Engine powering implementation of \c ECLResultData interface -class Opm::ECLResultData::Impl +// ====================================================================== +// Class (Anonymous)::ECLImpl::GridIDCache +// ====================================================================== + +ECLImpl::GridIDCache::GridIDCache(const ecl_file_view_type* block) + : block_(block) +{} + +int ECLImpl::GridIDCache::getGridID(const std::string& gridName) const +{ + if (gridName.empty()) { + return ECL_GRID_MAINGRID_LGR_NR; + } + + { + auto i = this->cache_.find(gridName); + if (i != std::end(this->cache_)) { + return i->second; + } + } + + const auto nLGR = ecl_file_view_get_num_named_kw + (this->block_, LGR_KW); + + auto lgrID = 0 * nLGR; + for (; lgrID < nLGR; ++lgrID) + { + const auto* kw = ecl_file_view_iget_named_kw + (this->block_, LGR_KW, lgrID); + + if ((getKeywordElementType(kw) != ECL_CHAR_TYPE) || + (ecl_kw_get_size(kw) != 1)) + { + // Huh !?! + continue; + } + + const auto kwname = GetKeywordData + ::as(kw, std::true_type()); + + if (kwname[0] == gridName) { + break; + } + } + + if (lgrID == nLGR) { + // No such LGR in block. Somewhat surprising. + return -1; + } + + return this->cache_[gridName] = + ECL_GRID_MAINGRID_LGR_NR + 1 + lgrID; +} + +// ====================================================================== +// Class (Anonymous)::ECLImpl::InitFileSections +// ====================================================================== + +ECLImpl::InitFileSections::InitFileSections(const ecl_file_type* init) + // Note: ecl_file_get_global_view() does not modify input arg + : init_(ecl_file_get_global_view(const_cast(init))) +{ + const auto* endLGR_kw = "LGRSGONE"; + const auto nEndLGR = + ecl_file_view_get_num_named_kw(this->init_, endLGR_kw); + + if (nEndLGR == 0) { + // No LGRs in model. INIT file consists of global section only, + // meaning that the only available section is equal to the global + // view (i.e., this->init_). + this->sect_.push_back(Section{ this->init_ }); + } + else { + const auto* start_kw = INTEHEAD_KW; + const auto* end_kw = endLGR_kw; + + this->sect_.reserve(2 * nEndLGR); + + for (auto sectID = 0*nEndLGR; sectID < nEndLGR; ++sectID) { + // Note: Start keyword occurrence lags one behind section ID + // when creating sections *between* LGRSGONE keywords. + const auto start_kw_occurrence = (sectID > 0) + ? sectID - 1 : 0*sectID; + + // Main section 'sectID': [ start_kw, LGRSGONE ] + const auto* sect = + ecl_file_view_add_blockview2(this->init_, start_kw, + end_kw, start_kw_occurrence); + + if (sect == nullptr) { + continue; + } + + const auto firstkw = Details::firstBlockKeyword(sect); + const auto occurrence = 0; + + // Main grid sub-section of 'sectID': [ start_kw, LGR ] + const auto* main_grid_sect = + ecl_file_view_add_blockview2(sect, firstkw.c_str(), + LGR_KW, occurrence); + + // LGR sub-section of 'sectID': [ LGR, LGRSGONE ] + const auto* lgr_sect = + ecl_file_view_add_blockview2(sect, LGR_KW, + end_kw, occurrence); + + // Note: main_grid_sect or lgr_sect *may* (in rare cases) be + // nullptr, but we'll deal with that in the calling context. + this->sect_.push_back(Section { main_grid_sect }); + this->sect_.push_back(Section { lgr_sect }); + + // start_kw == end_kw for all but first section. + start_kw = end_kw; + } + } +} + +// ====================================================================== +// Class Opm::ECLRestartData::Impl +// ====================================================================== + +/// Engine powering implementation of \c ECLRestartData interface +class Opm::ECLRestartData::Impl { public: + using Path = boost::filesystem::path; + /// Constructor /// - /// \param[in] rset Filesystem element or casename prefix representing + /// \param[in] rstrt Filesystem element or casename prefix representing /// an ECL result-set. - Impl(Path rset); + Impl(Path rstrt); + + /// Constructor + /// + /// \param[in] rstrt ECL restart result set + Impl(std::shared_ptr rstrt); /// Copy constructor. /// @@ -516,27 +826,6 @@ public: /// instance. Underlying result-set accessor is null upon return. Impl(Impl&& rhs); - /// Access the underlying ERT representation of the result-set. - /// - /// This is a hole in the interface that exists to be able to access - /// ERT's internal data structures for non-neighbouring connections - /// (i.e., faults and/or connections between main grid and LGRs). See - /// function ecl_nnc_export() in the ERT. - /// - /// 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. /// /// This is needed when working with dynamic restart data. @@ -560,7 +849,7 @@ public: /// \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& gridName) const; /// Retrieve current result-set view's data values for particular named /// result vector in particular enumerated grid. @@ -588,7 +877,7 @@ public: template std::vector keywordData(const std::string& vector, - const int gridID) const; + const std::string& gridName) const; private: /// RAII class to select a sub-block pertaining to a particular grid @@ -642,7 +931,7 @@ private: /// Saved original active view from host (prior to restricting view /// to single grid ID). - ecl_file_view_type* save_; + const ecl_file_view_type* save_; }; /// Casename prefix. Mostly to implement copy ctor. @@ -655,8 +944,11 @@ private: /// of main grid's section within a view. std::string firstKeyword_; + /// Map LGR names to integral grid IDs. + std::unique_ptr gridIDCache_; + /// Current active result-set view. - mutable ecl_file_view_type* activeBlock_{ nullptr }; + mutable const ecl_file_view_type* activeBlock_{ nullptr }; /// Support for passing \code *this \endcode to ERT functions that /// require an \c ecl_file_type, particularly the function that selects @@ -673,45 +965,42 @@ private: /// Retrieve result-set keyword that identifies beginning of main grid's /// result vectors. const std::string& mainGridStart() const; + + int gridID(const std::string& gridName) const; }; -Opm::ECLResultData::Impl::Impl(Path prefix) +Opm::ECLRestartData::Impl::Impl(Path prefix) : prefix_ (std::move(prefix)) - , result_ (openResultSet(deriveResultPath(prefix_))) + , result_ (openResultSet(deriveRestartPath(prefix_))) , firstKeyword_(firstFileKeyword(result_.get())) {} -Opm::ECLResultData::Impl::Impl(const Impl& rhs) +Opm::ECLRestartData::Impl::Impl(std::shared_ptr rstrt) + : prefix_ (ecl_file_get_src_file(rstrt.get())) + , result_ (std::move(rstrt)) + , firstKeyword_(firstFileKeyword(result_.get())) +{} + +Opm::ECLRestartData::Impl::Impl(const Impl& rhs) : prefix_ (rhs.prefix_) - , result_ (openResultSet(deriveResultPath(prefix_))) + , result_ (openResultSet(deriveRestartPath(prefix_))) , firstKeyword_(firstFileKeyword(result_.get())) {} -Opm::ECLResultData::Impl::Impl(Impl&& rhs) +Opm::ECLRestartData::Impl::Impl(Impl&& rhs) : prefix_ (std::move(rhs.prefix_)) , result_ (std::move(rhs.result_)) , firstKeyword_(std::move(rhs.firstKeyword_)) {} -const ecl_file_type* -Opm::ECLResultData::Impl::getRawFilePtr() const -{ - return this->result_.get(); -} - -bool Opm::ECLResultData::Impl::selectGlobalView() -{ - this->activeBlock_ = ecl_file_get_global_view(*this); - - return this->activeBlock_ != nullptr; -} - -bool Opm::ECLResultData::Impl::selectReportStep(const int step) +bool Opm::ECLRestartData::Impl::selectReportStep(const int step) { if (! ecl_file_has_report_step(*this, step)) { return false; } + this->gridIDCache_.reset(); + if (auto* globView = ecl_file_get_global_view(*this)) { // Ignore sequence numbers, dates, and simulation time. const auto seqnum = -1; @@ -722,16 +1011,26 @@ bool Opm::ECLResultData::Impl::selectReportStep(const int step) ecl_file_view_add_restart_view(globView, seqnum, step, dates, simdays); - return this->activeBlock_ != nullptr; + if (this->activeBlock_ != nullptr) { + this->gridIDCache_ + .reset(new ECLImpl::GridIDCache(this->activeBlock_)); + + return true; + } } return false; } -bool Opm::ECLResultData::Impl:: -haveKeywordData(const std::string& vector, const int gridID) const +bool Opm::ECLRestartData::Impl:: +haveKeywordData(const std::string& vector, + const std::string& gridName) const { - assert ((gridID >= 0) && "Grid IDs must be non-negative"); + const auto gridID = this->gridIDCache_->getGridID(gridName); + + if (gridID < 0) { + return false; + } // Note: Non-trivial dtor. Compiler can't ignore object. const auto block = Restrict{ *this, gridID }; @@ -746,18 +1045,22 @@ namespace Opm { template std::vector - ECLResultData::Impl::keywordData(const std::string& vector, - const int gridID) const + ECLRestartData::Impl::keywordData(const std::string& vector, + const std::string& gridName) const { - if (! this->haveKeywordData(vector, gridID)) { + if (! this->haveKeywordData(vector, gridName)) { std::ostringstream os; - os << "Cannot Access Non-Existent Keyword Data Pair (" - << vector << ", " << gridID << ')'; + os << "RESTART: Cannot Access Non-Existent Keyword Data Pair (" + << vector << ", " + << (gridName.empty() ? "Main Grid" : gridName) + << ')'; throw std::invalid_argument(os.str()); } + const auto gridID = this->gridIDCache_->getGridID(gridName); + // Note: Non-trivial dtor. Compiler can't ignore object. const auto block = Restrict{ *this, gridID }; @@ -770,105 +1073,483 @@ namespace Opm { assert ((kw != nullptr) && "Logic Error In Data Availability Check"); - // Whether or not caller requests a vector. - const auto makeStringVector = - typename std::is_same::type{}; - - switch (getKeywordElementType(kw)) { - case ECL_CHAR_TYPE: - return ECLImpl::GetKeywordData:: - as(kw, makeStringVector); - - case ECL_INT_TYPE: - return ECLImpl::GetKeywordData:: - as(kw, makeStringVector); - - case ECL_FLOAT_TYPE: - return ECLImpl::GetKeywordData:: - as(kw, makeStringVector); - - case ECL_DOUBLE_TYPE: - return ECLImpl::GetKeywordData:: - as(kw, makeStringVector); - - default: - // No operator exists for this type. Return empty. - return {}; - } + return ECLImpl::getKeywordData(kw); } } // namespace Opm -Opm::ECLResultData::Impl::operator ecl_file_type*() const +Opm::ECLRestartData::Impl::operator ecl_file_type*() const { return this->result_.get(); } -Opm::ECLResultData::Impl::operator const ecl_file_view_type*() const +Opm::ECLRestartData::Impl::operator const ecl_file_view_type*() const { return this->activeBlock_; } const std::string& -Opm::ECLResultData::Impl::mainGridStart() const +Opm::ECLRestartData::Impl::mainGridStart() const { return this->firstKeyword_; } +int +Opm::ECLRestartData::Impl::gridID(const std::string& gridName) const +{ + return this->gridIDCache_->getGridID(gridName); +} + // ====================================================================== -// Implementation of class Opm::ECLResultData Below Separator +// Class Opm::ECLInitFileData::Impl // ====================================================================== -Opm::ECLResultData::ECLResultData(const Path& prefix) - : pImpl_(new Impl(prefix)) +class Opm::ECLInitFileData::Impl +{ +public: + using Path = boost::filesystem::path; + + /// Constructor. + /// + /// Construct from filename. Owning semantics. + /// + /// \param[in] casePrefix Name or prefix of ECL result data. + Impl(Path initFile); + + /// Constructor. + /// + /// Construct from dataset already input through other means. + /// + /// Non-owning semantics. + Impl(std::shared_ptr initFile); + + /// Copy constructor. + /// + /// \param[in] rhs Object from which to construct new \c Impl instance. + Impl(const Impl& rhs); + + /// Move constructor. + /// + /// \param[in,out] rhs Object from which to constructo new \c Impl + /// instance. Underlying result-set accessor is null upon return. + Impl(Impl&& rhs); + + const ecl_file_type* getRawFilePtr() const; + + /// 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. + /// + /// \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& gridName) 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 + /// \endcode. Similarly, we cannot access underlying character data + /// through elements of an arithmetic type (e.g., \code + /// std::vector \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. + /// + /// \return Keyword data values. Empty if type conversion fails. + template + std::vector + keywordData(const std::string& vector, + const std::string& gridName) const; + +private: + using SectionID = + ECLImpl::InitFileSections::SectionID; + + /// Result of searching for a particular pairing of (vector,gridName) in + /// INIT. + struct LookupResult + { + /// In what INIT file section the KW was located (-1 if not found). + SectionID sectID; + + /// Local ID, within sectID, of grid section data section that hosts + /// the 'vector' (-1 if not found). + int gridSectID; + }; + + /// Pairing of kw vector and grid name. + struct KWKey { + /// Keyword/result set vector + std::string vector; + + /// ID of grid for which to look up 'vector'. + std::string gridName; + }; + + /// Comparator (std::set<> and std::map<>) for KWKeys. + struct CompareKWKey { + bool operator()(const KWKey& k1, const KWKey& k2) const + { + return std::tie(k1.vector, k1.gridName) + < std::tie(k2.vector, k2.gridName); + } + }; + + /// Negative look-up cache. + using MissingKW = std::set; + + /// Keyword-to-section cache to accelerate repeated look-up. + using KWSection = std::map; + + /// Casename prefix. Mostly to implement copy ctor. + const Path prefix_; + + /// Raw result set. + ECLImpl::FilePtr initFile_; + + /// Sections of the INIT result set. + ECLImpl::InitFileSections sections_; + + mutable const ecl_file_view_type* activeBlock_{ nullptr }; + + /// Negative look-up cache for haveKeywordData() queries. + mutable MissingKW missing_kw_; + + /// Keyword-to-section map for successful haveKeywordData() and + /// keywordData() queries. Accelerates repeated look-up queries. + mutable KWSection kw_section_; + + operator const ecl_file_view_type*() const; + + LookupResult + lookup(const std::string& vector, + const std::string& gridName) const; + + LookupResult lookupMainGrid(const KWKey& key) const; + + LookupResult lookupLGR(const KWKey& key) const; + + void setActiveBlock(const SectionID sect) const; + + const ECLImpl::InitFileSections::Section& + getSection(const SectionID sect) const; +}; + +Opm::ECLInitFileData::Impl::Impl(Path initFile) + : prefix_ (initFile.stem()) + , initFile_(openResultSet(deriveInitPath(std::move(initFile)))) + , sections_(initFile_.get()) {} -Opm::ECLResultData::ECLResultData(const ECLResultData& rhs) +Opm::ECLInitFileData::Impl::Impl(std::shared_ptr initFile) + : prefix_ (Path(ecl_file_get_src_file(initFile.get())).stem()) + , initFile_(std::move(initFile)) + , sections_(initFile_.get()) +{} + +Opm::ECLInitFileData::Impl::Impl(const Impl& rhs) + : prefix_ (rhs.prefix_) + , initFile_(rhs.initFile_) + , sections_(initFile_.get()) +{} + +Opm::ECLInitFileData::Impl::Impl(Impl&& rhs) + : prefix_ (std::move(rhs.prefix_)) + , initFile_(std::move(rhs.initFile_)) + , sections_(std::move(rhs.sections_)) +{} + +const ecl_file_type* +Opm::ECLInitFileData::Impl::getRawFilePtr() const +{ + return this->initFile_.get(); +} + +bool +Opm::ECLInitFileData::Impl:: +haveKeywordData(const std::string& vector, + const std::string& gridName) const +{ + const auto kwloc = this->lookup(vector, gridName); + + return (kwloc.sectID < this->sections_.numSections()) + && (kwloc.gridSectID >= 0); +} + +namespace Opm { + + template + std::vector + ECLInitFileData::Impl:: + keywordData(const std::string& vector, + const std::string& gridName) const + { + if (! this->haveKeywordData(vector, gridName)) { + std::ostringstream os; + + os << "INIT: Cannot Access Non-Existent Keyword Data Pair (" + << vector << ", " + << (gridName.empty() ? "Main Grid" : gridName) + << ')'; + + throw std::invalid_argument(os.str()); + }; + + const auto kwloc = this->lookup(vector, gridName); + + this->setActiveBlock(kwloc.sectID); + + if (! gridName.empty()) { + // We're cons + this->activeBlock_ = + ecl_file_view_add_blockview(this->activeBlock_, LGR_KW, + kwloc.gridSectID); + } + + const auto occurrence = 0; + + const auto* kw = + ecl_file_view_iget_named_kw(*this, vector.c_str(), + occurrence); + + assert ((kw != nullptr) && + "Logic Error In Data Availability Check"); + + return ECLImpl::getKeywordData(kw); + } + +} + +Opm::ECLInitFileData::Impl::operator const ecl_file_view_type*() const +{ + return this->activeBlock_; +} + +Opm::ECLInitFileData::Impl::LookupResult +Opm::ECLInitFileData::Impl:: +lookup(const std::string& vector, const std::string& gridName) const +{ + const auto key = KWKey{ vector, gridName }; + + { + auto m = this->missing_kw_.find(key); + + if (m != std::end(this->missing_kw_)) { + // 'vector' known to be mssing for 'gridName'. Report as such. + return { SectionID(-1), -1 }; + } + } + + { + auto i = this->kw_section_.find(key); + + if (i != std::end(this->kw_section_)) { + // 'vector' previously located for 'gridName'. Return it. + return i->second; + } + } + + if (gridName.empty()) { + return this->lookupMainGrid(key); + } + + return this->lookupLGR(key); +} + +Opm::ECLInitFileData::Impl::LookupResult +Opm::ECLInitFileData::Impl::lookupMainGrid(const KWKey& key) const +{ + assert (key.gridName.empty() && "Logic Error"); + + const auto* kwheader = key.vector.c_str(); + + // Main grid sections are even numbered. + for (auto nSect = this->sections_.numSections(), sectID = 0*nSect; + sectID < nSect; sectID += 2) + { + const auto& s = this->getSection(sectID); + + // Skip empty sections + if (s.block == nullptr) { continue; } + + const auto count = + ecl_file_view_get_num_named_kw(s.block, kwheader); + + if (count > 0) { + // Result-set 'vector' present in main grid. Record and return. + + // Assume that keyword does not occur multiple times in main + // grid section. + const auto gridSectID = 0; + + return this->kw_section_[key] = + LookupResult { sectID, gridSectID }; + } + } + + // No result-set 'vector' in main grid. Record as missing and return + // "not found". + this->missing_kw_.insert(key); + + return { SectionID(-1), -1 }; +} + +Opm::ECLInitFileData::Impl::LookupResult +Opm::ECLInitFileData::Impl::lookupLGR(const KWKey& key) const +{ + assert ((! key.gridName.empty()) && "Logic Error"); + + const auto gridID = paddedGridName(key.gridName); + const auto* kwheader = key.vector.c_str(); + const auto* gridName = gridID.c_str(); + + // LGR sections are odd numbered. + for (auto nSect = this->sections_.numSections(), sectID = 0*nSect + 1; + sectID < nSect; sectID += 2) + { + const auto& s = this->getSection(sectID); + + // Skip empty sections. + if (s.block == nullptr) { continue; } + + const auto nLGR = + ecl_file_view_get_num_named_kw(s.block, LGR_KW); + + assert ((nLGR > 0) && "Logic Error"); + + for (auto lgrID = 0*nLGR; lgrID < nLGR; ++lgrID) { + const auto* kw = + ecl_file_view_iget_named_kw(s.block, LGR_KW, lgrID); + + if (! ecl_kw_data_equal(kw, gridName)) { + // This is not the grid you're looking for. + continue; + } + + // This LGR section pertains to key.gridName. Look for + // key.vector between this occurrence of LGR_KW and the next. + // Continue searching if we don't find that vector however, + // because there may be multiple relevant LGR instances for + // key.gridName. + + const auto* lgrSect = + ecl_file_view_add_blockview(s.block, LGR_KW, lgrID); + + const auto kwCount = + ecl_file_view_get_num_named_kw(lgrSect, kwheader); + + if (kwCount > 0) { + // We found the key.vector in this LGR data section. Record + // position and return it to caller. + + return this->kw_section_[key] = + LookupResult { SectionID(sectID), lgrID }; + } + + // Did not find 'key.vector' in this LGR data section. Continue + // searching nonetheless because there may be more LGR sections + // that pertain to 'key.gridName'. + } + } + + // key.vector not found for key.gridName in any of the LGR sections. + // Record as missing and return "not found". + + this->missing_kw_.insert(key); + + return { SectionID(-1), -1 }; +} + +void +Opm::ECLInitFileData::Impl::setActiveBlock(const SectionID sect) const +{ + this->activeBlock_ = this->getSection(sect).block; +} + +const ECLImpl::InitFileSections::Section& +Opm::ECLInitFileData::Impl::getSection(const SectionID sect) const +{ + assert (sect < this->sections_.numSections()); + + return this->sections_[sect]; +} + +// ###################################################################### +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = +// Public Interfaces Follow +// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// ###################################################################### + +// ====================================================================== +// Implementation of class Opm::ECLRestartData Below Separator +// ====================================================================== + +Opm::ECLRestartData::ECLRestartData(boost::filesystem::path rstrt) + : pImpl_(new Impl(std::move(rstrt))) +{} + +Opm::ECLRestartData::ECLRestartData(std::shared_ptr rstrt) + : pImpl_(new Impl(std::move(rstrt))) +{} + +Opm::ECLRestartData::ECLRestartData(const ECLRestartData& rhs) : pImpl_(new Impl(*rhs.pImpl_)) {} -Opm::ECLResultData::ECLResultData(ECLResultData&& rhs) +Opm::ECLRestartData::ECLRestartData(ECLRestartData&& rhs) : pImpl_(std::move(rhs.pImpl_)) {} -Opm::ECLResultData& -Opm::ECLResultData::operator=(const ECLResultData& rhs) +Opm::ECLRestartData& +Opm::ECLRestartData::operator=(const ECLRestartData& rhs) { this->pImpl_.reset(new Impl(*rhs.pImpl_)); return *this; } -Opm::ECLResultData& -Opm::ECLResultData::operator=(ECLResultData&& rhs) +Opm::ECLRestartData& +Opm::ECLRestartData::operator=(ECLRestartData&& rhs) { this->pImpl_ = std::move(rhs.pImpl_); return *this; } -Opm::ECLResultData::~ECLResultData() +Opm::ECLRestartData::~ECLRestartData() {} -const ecl_file_type* -Opm::ECLResultData::getRawFilePtr() const +bool Opm::ECLRestartData::selectReportStep(const int step) const { - return this->pImpl_->getRawFilePtr(); -} + // selectReportStep() const is a bit of a lie. pImpl_ is a const + // pointer to modifiable Impl. -bool Opm::ECLResultData::selectGlobalView() -{ - return this->pImpl_->selectGlobalView(); -} - -bool Opm::ECLResultData::selectReportStep(const int step) -{ return this->pImpl_->selectReportStep(step); } bool -Opm::ECLResultData:: -haveKeywordData(const std::string& vector, const int gridID) const +Opm::ECLRestartData:: +haveKeywordData(const std::string& vector, + const std::string& gridID) const { return this->pImpl_->haveKeywordData(vector, gridID); } @@ -877,23 +1558,101 @@ namespace Opm { template std::vector - ECLResultData::keywordData(const std::string& vector, - const int gridID) const + ECLRestartData::keywordData(const std::string& vector, + const std::string& gridID) const { return this->pImpl_->template keywordData(vector, gridID); } // Explicit instantiations for those types we care about. template std::vector - ECLResultData::keywordData(const std::string& vector, - const int gridID) const; + ECLRestartData::keywordData(const std::string& vector, + const std::string& gridID) const; template std::vector - ECLResultData::keywordData(const std::string& vector, - const int gridID) const; + ECLRestartData::keywordData(const std::string& vector, + const std::string& gridID) const; template std::vector - ECLResultData::keywordData(const std::string& vector, - const int gridID) const; + ECLRestartData::keywordData(const std::string& vector, + const std::string& gridID) const; + +} // namespace Opm::ECL + +// ====================================================================== +// Implementation of class Opm::ECLInitFileData Below Separator +// ====================================================================== + +Opm::ECLInitFileData::ECLInitFileData(boost::filesystem::path init) + : pImpl_(new Impl(std::move(init))) +{} + +Opm::ECLInitFileData::ECLInitFileData(std::shared_ptr init) + : pImpl_(new Impl(std::move(init))) +{} + +Opm::ECLInitFileData::ECLInitFileData(const ECLInitFileData& rhs) + : pImpl_(new Impl(*rhs.pImpl_)) +{} + +Opm::ECLInitFileData::ECLInitFileData(ECLInitFileData&& rhs) + : pImpl_(std::move(rhs.pImpl_)) +{} + +Opm::ECLInitFileData& +Opm::ECLInitFileData::operator=(const ECLInitFileData& rhs) +{ + this->pImpl_.reset(new Impl(*rhs.pImpl_)); + + return *this; +} + +Opm::ECLInitFileData& +Opm::ECLInitFileData::operator=(ECLInitFileData&& rhs) +{ + this->pImpl_ = std::move(rhs.pImpl_); + + return *this; +} + +Opm::ECLInitFileData::~ECLInitFileData() +{} + +bool +Opm::ECLInitFileData:: +haveKeywordData(const std::string& vector, + const std::string& gridID) const +{ + return this->pImpl_->haveKeywordData(vector, gridID); +} + +const ecl_file_type* +Opm::ECLInitFileData::getRawFilePtr() const +{ + return this->pImpl_->getRawFilePtr(); +} + +namespace Opm { + + template + std::vector + ECLInitFileData::keywordData(const std::string& vector, + const std::string& gridID) const + { + return this->pImpl_->template keywordData(vector, gridID); + } + + // Explicit instantiations for those types we care about. + template std::vector + ECLInitFileData::keywordData(const std::string& vector, + const std::string& gridID) const; + + template std::vector + ECLInitFileData::keywordData(const std::string& vector, + const std::string& gridID) const; + + template std::vector + ECLInitFileData::keywordData(const std::string& vector, + const std::string& gridID) const; } // namespace Opm::ECL diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.hpp index a804e4baf8..c48b6a4645 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.hpp @@ -1,6 +1,5 @@ /* - 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). @@ -18,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef OPM_ECLRESTARTDATA_HEADER_INCLUDED -#define OPM_ECLRESTARTDATA_HEADER_INCLUDED +#ifndef OPM_ECLRESULTDATA_HEADER_INCLUDED +#define OPM_ECLRESULTDATA_HEADER_INCLUDED #include #include @@ -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. - explicit ECLResultData(const Path& casePrefix); + /// Owning semantics. + /// + /// \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 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(); - - /// 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(); + ~ECLRestartData(); /// 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 std::vector 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 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 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 \endcode. Similarly, we cannot + /// access underlying character data through elements of an + /// arithmetic type (e.g., \code std::vector \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 + std::vector + 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 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 diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.cpp index 4a3b77b1c1..4f29ba1ae6 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.cpp +++ b/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). diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.hpp index ee11974fe0..87fe95d75b 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.hpp +++ b/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). diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp index 20b611a4a0..f1c17fbc52 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp +++ b/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::solution(const ECLResultData& restart, - const int num_grids) const + ECLWellSolution::solution(const ECLRestartData& restart, + const std::vector& 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 all_wd{}; - for (int grid_index = 0; grid_index < num_grids; ++grid_index) { - std::vector wd = readWellData(restart, grid_index); + for (const auto& gridName : grids) { + std::vector 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::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 - // where the correct restart block using the ert block mechanisms. + // Check if result set provides complete set of well solution data. + 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(INTEHEAD_KW, grid_index)); + INTEHEAD ih(restart.keywordData(INTEHEAD_KW, gridName)); if (ih.nwell == 0) { return {}; } const double qr_unit = resRateUnit(ih.unit); - // Read necessary keywords. - auto zwel = restart.keywordData(ZWEL_KW, grid_index); - auto iwel = restart.keywordData (IWEL_KW, grid_index); - auto xwel = restart.keywordData ("XWEL" , grid_index); - auto icon = restart.keywordData (ICON_KW, grid_index); - auto xcon = restart.keywordData ("XCON" , grid_index); + // Load well topology and flow rates. + auto zwel = restart.keywordData(ZWEL_KW, gridName); + auto iwel = restart.keywordData (IWEL_KW, gridName); + auto xwel = restart.keywordData ("XWEL" , gridName); + auto icon = restart.keywordData (ICON_KW, gridName); + auto xcon = restart.keywordData ("XCON" , gridName); // Create well data. std::vector 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 { diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp index e6ab682b2a..94e24c65ec 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp +++ b/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 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 solution(const ECLResultData& restart, - const int num_grids) const; + std::vector solution(const ECLRestartData& restart, + const std::vector& grids) const; private: // Data members. @@ -67,8 +67,8 @@ namespace Opm bool disallow_crossflow_; // Methods. - std::vector readWellData(const ECLResultData& restart, - const int grid_index) const; + std::vector readWellData(const ECLRestartData& restart, + const std::string& gridName) const; }; diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runAcceptanceTest.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runAcceptanceTest.cpp index 3b799a11fb..bb8108d72d 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runAcceptanceTest.cpp +++ b/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) { diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runLinearisedCellDataTest.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runLinearisedCellDataTest.cpp index 89500d9b8e..1811f8a3ed 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runLinearisedCellDataTest.cpp +++ b/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 - testQuantities(const ::Opm::parameter::ParameterGroup& param) + testQuantities(const ::Opm::ParameterGroup& param) { return StringUtils::VectorValue:: get(param.get("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 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& steps) { @@ -510,7 +511,7 @@ namespace { auto E = std::array{}; 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(ECLquant) + graph.rawLinearisedCellData(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); diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runTransTest.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runTransTest.cpp index e705cfede4..766a11048c 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runTransTest.cpp +++ b/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 - 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& trans) { const auto Tref = loadReference(param); diff --git a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.cpp b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.cpp index 3ad859a0ef..dca1df6d21 100644 --- a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.cpp +++ b/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& pv) { - const auto& ftof = injector_solution.fd.timeOfFlight(); - const auto& rtof = producer_solution.fd.timeOfFlight(); - if (pv.size() != ftof.size() || pv.size() != rtof.size()) { + return flowCapacityStorageCapacityCurve(injector_solution.fd.timeOfFlight(), + producer_solution.fd.timeOfFlight(), + 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& injector_tof, + const std::vector& producer_tof, + const std::vector& 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 D2; std::vector 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)}; diff --git a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.hpp b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.hpp index 0f3642bec1..22023ab708 100644 --- a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.hpp +++ b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.hpp @@ -50,6 +50,14 @@ namespace FlowDiagnostics const Toolbox::Reverse& producer_solution, const std::vector& 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& injector_tof, + const std::vector& producer_tof, + const std::vector& pore_volume); + /// The Lorenz coefficient from the F-Phi curve. /// /// The Lorenz coefficient is a measure of heterogeneity. It diff --git a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/Solution.cpp b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/Solution.cpp index ee88861d39..d31e65bc7e 100644 --- a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/Solution.cpp +++ b/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). diff --git a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/TracerTofSolver.cpp b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/TracerTofSolver.cpp index 064f8c61e7..dc1994df86 100644 --- a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/TracerTofSolver.cpp +++ b/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). diff --git a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/utility/graph/AssembledConnectionsIteration.hpp b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/utility/graph/AssembledConnectionsIteration.hpp index b949242505..d91baeaeff 100644 --- a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/utility/graph/AssembledConnectionsIteration.hpp +++ b/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). diff --git a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_derivedquantities.cpp b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_derivedquantities.cpp index bab38628a1..6c14694949 100644 --- a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_derivedquantities.cpp +++ b/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 pv { 1.0, 2.0, 1.0 }; + std::vector ftof { 0.0, 2.0, 1.0 }; + std::vector 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()