From 8596c3f79c0a68c0affe918cd067df7b96a34efd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jacob=20St=C3=B8ren?= Date: Thu, 9 Feb 2017 14:14:02 +0100 Subject: [PATCH] #1206 Updated the opm-flowdiag libraries --- ResInsightVersion.cmake | 6 +- .../AcceptanceTests.cmake | 58 ++ .../CMakeLists.txt | 8 + .../CMakeLists_files.cmake | 8 + .../examples/exampleSetup.hpp | 155 ++--- .../opm/utility/ECLFluxCalc.cpp | 64 ++ .../opm/utility/ECLFluxCalc.hpp | 61 ++ .../opm/utility/ECLGraph.cpp | 505 +++++++++++++-- .../opm/utility/ECLGraph.hpp | 69 +- .../opm/utility/ECLUnitHandling.cpp | 209 +++++++ .../opm/utility/ECLUnitHandling.hpp | 45 ++ .../opm/utility/ECLWellSolution.cpp | 11 +- .../tests/runAcceptanceTest.cpp | 550 ++++++++++++++++ .../tests/runLinearisedCellDataTest.cpp | 587 ++++++++++++++++++ .../tests/runTransTest.cpp | 229 +++++++ .../tests/test_eclunithandling.cpp | 235 +++++++ .../opm/parser/eclipse/Units/Units.hpp | 445 +++++++------ .../opm/flowdiagnostics/TracerTofSolver.cpp | 10 +- 18 files changed, 2945 insertions(+), 310 deletions(-) create mode 100644 ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/AcceptanceTests.cmake create mode 100644 ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp create mode 100644 ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp create mode 100644 ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.cpp create mode 100644 ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.hpp create mode 100644 ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runAcceptanceTest.cpp create mode 100644 ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runLinearisedCellDataTest.cpp create mode 100644 ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runTransTest.cpp create mode 100644 ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclunithandling.cpp diff --git a/ResInsightVersion.cmake b/ResInsightVersion.cmake index 2008251a16..aead5be760 100644 --- a/ResInsightVersion.cmake +++ b/ResInsightVersion.cmake @@ -13,8 +13,10 @@ set(ERT_GITHUB_SHA "e2a5a9cc20705537d07822958d925e092a323367") set(OPM_COMMON_GITHUB_SHA "1216bc052542f24ec6fcfbe1947d52e6300ff754") set(OPM_PARSER_GITHUB_SHA "a3496df501a4369fd827fbf0ff893d03deff425f") -set(OPM_FLOWDIAGNOSTICS_SHA "80190c28ae0fd4a82cfe88c827559029982db83b") -set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "01ecb8fc22cb70d883edaad398bffc49878859c3") +set(OPM_FLOWDIAGNOSTICS_SHA "02503abf0043bb08062e358e460a0c8231b6225c") +set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "09057e5767e492fa44c5e9ae78e9e1fcc9694b6c") + +# sha for Units.hpp 9a679071dd0066236154852c39a9e0b3c3ac4873 set(STRPRODUCTVER ${RESINSIGHT_MAJOR_VERSION}.${RESINSIGHT_MINOR_VERSION}.${RESINSIGHT_INCREMENT_VERSION}) diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/AcceptanceTests.cmake b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/AcceptanceTests.cmake new file mode 100644 index 0000000000..debbfb1856 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/AcceptanceTests.cmake @@ -0,0 +1,58 @@ +# Set absolute tolerance to be used for testing +Set (abs_tol 5.0e-8) +Set (rel_tol 1.0e-13) + +# Input: +# - casename: with or without extension +# +Macro (add_acceptance_test casename) + + String (REGEX REPLACE "\\.[^.]*$" "" basename "${casename}") + + Add_Test (NAME ToF_accept_${casename}_all_steps + COMMAND runAcceptanceTest + "case=${OPM_DATA_ROOT}/flow_diagnostic_test/eclipse-simulation/${basename}" + "ref-dir=${OPM_DATA_ROOT}/flow_diagnostic_test/fd-ref-data/${basename}" + "atol=${abs_tol}" "rtol=${rel_tol}") + +EndMacro (add_acceptance_test) + +# Input +# - casename: with or without extension +# +Macro (add_trans_acceptance_test casename) + + String (REGEX REPLACE "\\.[^.]*$" "" basename "${casename}") + + Add_Test (NAME Trans_accept_${casename} + COMMAND runTransTest + "case=${OPM_DATA_ROOT}/flow_diagnostic_test/eclipse-simulation/${basename}" + "ref-dir=${OPM_DATA_ROOT}/flow_diagnostic_test/fd-ref-data/${basename}" + "atol=${abs_tol}" "rtol=${rel_tol}") + +EndMacro (add_trans_acceptance_test) + +# Input +# - casename: with or without extension +# - strings identifying which physical quantities to compare +Macro (add_celldata_acceptance_test casename) + + String (REGEX REPLACE "\\.[^.]*$" "" basename "${casename}") + + Add_Test (NAME CellData_accept_${casename} + COMMAND runLinearisedCellDataTest + "case=${OPM_DATA_ROOT}/flow_diagnostic_test/eclipse-simulation/${basename}" + "ref-dir=${OPM_DATA_ROOT}/flow_diagnostic_test/fd-ref-data/${basename}" + "quant=${ARGN}" "atol=${abs_tol}" "rtol=${rel_tol}") + +EndMacro (add_celldata_acceptance_test) + +If (NOT TARGET test-suite) + Add_Custom_Target (test-suite) +EndIf () + +# Acceptance tests + +Add_Acceptance_Test (SIMPLE_2PH_W_FAULT_LGR) +Add_Trans_Acceptance_Test (SIMPLE_2PH_W_FAULT_LGR) +Add_CellData_Acceptance_Test (SIMPLE_2PH_W_FAULT_LGR "pressure") diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists.txt b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists.txt index d6ca46957f..482918b86e 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists.txt +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists.txt @@ -59,6 +59,10 @@ endif() macro (dir_hook) endmacro (dir_hook) +# Look for the opm-data repository; if found the variable +# HAVE_OPM_DATA will be set to true. +include (Findopm-data) + # list of prerequisites for this particular project; this is in a # separate file (in cmake/Modules sub-directory) because it is shared # with the find module @@ -91,3 +95,7 @@ endmacro (install_hook) # all setup common to the OPM library modules is done here include (OpmLibMain) + +if (HAVE_OPM_DATA) + include (${CMAKE_CURRENT_SOURCE_DIR}/AcceptanceTests.cmake) +endif() 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 525d6de6b4..a75e534b84 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 @@ -21,22 +21,30 @@ # the library needs it. list (APPEND MAIN_SOURCE_FILES + opm/utility/ECLFluxCalc.cpp opm/utility/ECLGraph.cpp opm/utility/ECLResultData.cpp + opm/utility/ECLUnitHandling.cpp opm/utility/ECLWellSolution.cpp ) list (APPEND TEST_SOURCE_FILES + tests/test_eclunithandling.cpp ) list (APPEND EXAMPLE_SOURCE_FILES examples/computeLocalSolutions.cpp examples/computeToFandTracers.cpp examples/computeTracers.cpp + tests/runAcceptanceTest.cpp + tests/runLinearisedCellDataTest.cpp + tests/runTransTest.cpp ) list (APPEND PUBLIC_HEADER_FILES + opm/utility/ECLFluxCalc.hpp opm/utility/ECLGraph.hpp opm/utility/ECLResultData.hpp + opm/utility/ECLUnitHandling.hpp opm/utility/ECLWellSolution.hpp ) 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 389a4f2b4f..ab62d6ff88 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 @@ -29,6 +29,7 @@ #include #include +#include #include #include @@ -80,38 +81,31 @@ namespace example { } inline Opm::FlowDiagnostics::ConnectionValues - extractFluxField(const Opm::ECLGraph& G) + extractFluxField(const Opm::ECLGraph& G, const bool compute_fluxes) { using ConnVals = Opm::FlowDiagnostics::ConnectionValues; - - using NConn = ConnVals::NumConnections; - using NPhas = ConnVals::NumPhases; - - const auto nconn = NConn{G.numConnections()}; - const auto nphas = NPhas{3}; - - auto flux = ConnVals(nconn, nphas); + auto flux = ConnVals(ConnVals::NumConnections{G.numConnections()}, + ConnVals::NumPhases{3}); auto phas = ConnVals::PhaseID{0}; - for (const auto& p : { Opm::ECLGraph::PhaseIndex::Aqua , - Opm::ECLGraph::PhaseIndex::Liquid , - Opm::ECLGraph::PhaseIndex::Vapour }) + Opm::ECLFluxCalc calc(G); + + const auto phases = { Opm::ECLGraph::PhaseIndex::Aqua , + Opm::ECLGraph::PhaseIndex::Liquid , + Opm::ECLGraph::PhaseIndex::Vapour }; + + for (const auto& p : phases) { - const auto pflux = G.flux(p); - + const auto pflux = compute_fluxes ? calc.flux(p) : G.flux(p); if (! pflux.empty()) { - assert (pflux.size() == nconn.total); - + 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; } @@ -138,72 +132,60 @@ namespace example { return inflow; } - namespace Hack { - inline Opm::FlowDiagnostics::ConnectionValues - convert_flux_to_SI(Opm::FlowDiagnostics::ConnectionValues&& fl) + + + + struct FilePaths + { + FilePaths(const Opm::parameter::ParameterGroup& param) { - using Co = Opm::FlowDiagnostics::ConnectionValues::ConnID; - using Ph = Opm::FlowDiagnostics::ConnectionValues::PhaseID; - - const auto nconn = fl.numConnections(); - const auto nphas = fl.numPhases(); - - for (auto phas = Ph{0}; phas.id < nphas; ++phas.id) { - for (auto conn = Co{0}; conn.id < nconn; ++conn.id) { - fl(conn, phas) /= 86400; - } - } - - return fl; + const string casename = param.getDefault("case", "DEFAULT_CASE_NAME"); + grid = param.has("grid") ? param.get("grid") + : deriveFileName(casename, { ".EGRID", ".FEGRID", ".GRID", ".FGRID" }); + init = param.has("init") ? param.get("init") + : deriveFileName(casename, { ".INIT", ".FINIT" }); + restart = param.has("restart") ? param.get("restart") + : deriveFileName(casename, { ".UNRST", ".FUNRST" }); } - } - inline Opm::ECLGraph - initGraph(int argc, char* argv[]) + using path = boost::filesystem::path; + using string = std::string; + + path grid; + path init; + path restart; + }; + + + + + inline Opm::parameter::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); + return param; + } - // Obtain filenames for grid, init and restart files, as well as step number. - using boost::filesystem::path; - using std::string; - const string casename = param.getDefault("case", "DEFAULT_CASE_NAME"); - const path grid = param.has("grid") ? param.get("grid") - : deriveFileName(casename, { ".EGRID", ".FEGRID", ".GRID", ".FGRID" }); - const path init = param.has("init") ? param.get("init") - : deriveFileName(casename, { ".INIT", ".FINIT" }); - const path restart = param.has("restart") ? param.get("restart") - : deriveFileName(casename, { ".UNRST", ".FUNRST" }); - const int step = param.getDefault("step", 0); - // Read graph and fluxes, initialise the toolbox. - auto graph = Opm::ECLGraph::load(grid, init); - graph.assignFluxDataSource(restart); - if (! graph.selectReportStep(step)) { - std::ostringstream os; - - os << "Report Step " << step - << " is Not Available in Result Set '" - << grid.stem() << '\''; - - throw std::domain_error(os.str()); - } + 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; } - inline std::vector - initWellFluxes(const Opm::ECLGraph& G) - { - auto wsol = Opm::ECLWellSolution{}; - return wsol.solution(G.rawResultData(), G.numGrids()); - } + + inline Opm::FlowDiagnostics::Toolbox - initToolbox(const Opm::ECLGraph& G, const std::vector& well_fluxes) + initToolbox(const Opm::ECLGraph& G) { const auto connGraph = Opm::FlowDiagnostics:: ConnectivityGraph{ static_cast(G.numCells()), @@ -211,11 +193,7 @@ namespace example { // Create the Toolbox. auto tool = Opm::FlowDiagnostics::Toolbox{ connGraph }; - tool.assignPoreVolume(G.poreVolume()); - tool.assignConnectionFlux(Hack::convert_flux_to_SI(extractFluxField(G))); - - tool.assignInflowFlux(extractWellFlows(G, well_fluxes)); return tool; } @@ -226,15 +204,42 @@ namespace example { struct Setup { Setup(int argc, char** argv) - : graph(initGraph(argc, argv)) - , well_fluxes(initWellFluxes(graph)) - , toolbox(initToolbox(graph, well_fluxes)) + : param(initParam(argc, argv)) + , file_paths(param) + , 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)) { + std::ostringstream os; + os << "Report Step " << step + << " 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 { + return false; + } + } + + Opm::parameter::ParameterGroup param; + FilePaths file_paths; Opm::ECLGraph graph; std::vector well_fluxes; Opm::FlowDiagnostics::Toolbox toolbox; + bool compute_fluxes_ = false; }; 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 new file mode 100644 index 0000000000..b29636a6c8 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp @@ -0,0 +1,64 @@ +/* + Copyright 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 . +*/ + +#include + +namespace Opm +{ + + ECLFluxCalc::ECLFluxCalc(const ECLGraph& graph) + : graph_(graph) + , neighbours_(graph.neighbours()) + , transmissibility_(graph.transmissibility()) + { + } + + + + + + std::vector ECLFluxCalc::flux(const PhaseIndex /* phase */) const + { + // Obtain pressure vector. + std::vector pressure = graph_.linearisedCellData("PRESSURE", &ECLUnits::UnitSystem::pressure); + + // Compute fluxes per connection. + const int num_conn = transmissibility_.size(); + std::vector fluxvec(num_conn); + for (int conn = 0; conn < num_conn; ++conn) { + fluxvec[conn] = singleFlux(conn, pressure); + } + return fluxvec; + } + + + + + + double ECLFluxCalc::singleFlux(const int connection, + const std::vector& pressure) const + { + const int c1 = neighbours_[2*connection]; + const int c2 = neighbours_[2*connection + 1]; + const double t = transmissibility_[connection]; + return t * (pressure[c1] - pressure[c2]); + } + + +} // namespace Opm 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 new file mode 100644 index 0000000000..76a5409e93 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp @@ -0,0 +1,61 @@ +/* + Copyright 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 . +*/ + +#ifndef OPM_ECLFLUXCALC_HEADER_INCLUDED +#define OPM_ECLFLUXCALC_HEADER_INCLUDED + +#include +#include + +namespace Opm +{ + + /// Class for computing connection fluxes in the absence of flux output. + class ECLFluxCalc + { + public: + /// Construct from ECLGraph. + /// + /// \param[in] graph Connectivity data, as well as providing a means to read data from the restart file. + explicit ECLFluxCalc(const ECLGraph& graph); + + using PhaseIndex = ECLGraph::PhaseIndex; + + /// Retrive phase flux on all connections defined by \code + /// graph.neighbours() \endcode. + /// + /// \param[in] phase Canonical phase for which to retrive flux. + /// + /// \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; + + private: + double singleFlux(const int connection, + const std::vector& pressure) const; + + const ECLGraph& graph_; + std::vector neighbours_; + std::vector transmissibility_; + }; + +} // namespace Opm + +#endif // OPM_ECLFLUXCALC_HEADER_INCLUDED 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 e552b8991c..d06645c93d 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 @@ -24,6 +24,9 @@ #include #include +#include + +#include #include #include @@ -41,6 +44,7 @@ #include #include +#include #include #include @@ -93,6 +97,18 @@ namespace { std::array cartesianDimensions(const ecl_grid_type* G); + /// Access unit conventions pertaining to single grid in result set. + /// + /// \param[in] rset Result set. + /// + /// \param[in] grid_ID Numerical ID of grid. Non-negative. Zero + /// for the main grid and positive for LGRs. + /// + /// \return Unit system convention for \p grid_ID in result set. + auto getUnitSystem(const ::Opm::ECLResultData& rset, + const int grid_ID) + -> decltype(::Opm::ECLUnits::createUnitSystem(0)); + /// Retrieve global pore-volume vector from INIT source. /// /// Specialised tool needed to determine the active cells. @@ -101,7 +117,11 @@ namespace { /// /// \param[in] init ERT representation of INIT source. /// - /// \return Vector of pore-volumes for all global cells of \p G. + /// \param[in] grid_ID Numerical ID of grid. Non-negative. Zero + /// for the main grid and positive 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, @@ -137,6 +157,11 @@ namespace { const ::Opm::ECLResultData& init, const int gridID); + /// Retrieve non-negative numeric ID of grid instance. + /// + /// \return Constructor's \c gridID parameter. + int gridID() const; + /// Retrieve number of active cells in graph. std::size_t numCells() const; @@ -154,9 +179,18 @@ namespace { /// /// Corresponds to the \c PORV vector in the INIT file, possibly /// restricted to those active cells for which the pore-volume is - /// strictly positive. + /// strictly positive. SI unit conventions (rm^3). const std::vector& activePoreVolume() const; + /// Retrieve static (background) transmissibility values on all + /// connections defined by \code neighbours() \endcode. + /// + /// Specifically, \code transmissibility()[i] \endcode is the + /// transmissibility of the connection between cells \code + /// neighbours()[2*i + 0] \endcode and \code neighbours()[2*i + + /// 1] \endcode. + const std::vector& transmissibility() const; + /// Retrieve ID of active cell from global ID. int activeCell(const std::size_t globalCell) const; @@ -188,6 +222,11 @@ namespace { cellData(const ::Opm::ECLResultData& src, const std::string& vector) const; + template + std::vector + activeCellData(const ::Opm::ECLResultData& src, + const std::string& vector) const; + /// Retrieve values of result set vector for all Cartesian /// connections in grid. /// @@ -200,7 +239,8 @@ namespace { /// all of the grid's Cartesian connections. std::vector connectionData(const ::Opm::ECLResultData& src, - const std::string& vector) const; + const std::string& vector, + const double unit) const; private: /// Facility for deriving Cartesian neighbourship in a grid @@ -217,14 +257,18 @@ namespace { /// \param[in] G ERT Grid representation. /// /// \param[in] pvol Vector of pore-volumes on all global - /// cells of \p G. Typically obtained - /// through function getPVolVector(). + /// cells of \p G. Typically obtained through function + /// getPVolVector(). Numerical values assumed to be in + /// SI units (rm^3). CartesianCells(const ecl_grid_type* G, const std::vector& pvol); /// Retrive global cell indices of all active cells in grid. std::vector activeGlobal() const; + /// Retrieve pore-volume values for all active cells in grid. + /// + /// SI unit conventions (rm^3). const std::vector& activePoreVolume() const; /// Map input vector to all global cells. @@ -240,6 +284,30 @@ namespace { std::vector scatterToGlobal(const std::vector& x) const; + /// Restrict input vector to active grid cells. + /// + /// Selects subsets corresponding to active cells (i.e., + /// those cells for which \code pore_volume > 0 \endcode) if + /// input size is + /// + /// - All global cells + /// - Explicitly active cells (ACTNUM != 0) + /// + /// All other cases are returned unfiltered--i.e., as direct + /// copies of the input. + /// + /// \param[in] x Input vector, defined on the explicitly + /// active cells, all global cells or some + /// other subset (e.g., all non-neighbouring + /// connections). + /// + /// \return Input vector restricted to active cells if + /// subset known. Direct copy if \p x is defined on set + /// other than explicitly active or all global cells. + template + std::vector + gatherToActive(const std::vector& x) const; + /// Retrieve total number of cells in grid, including /// inactive ones. /// @@ -315,13 +383,8 @@ namespace { /// Whether or not a particular active cell is subdivided. std::vector is_divided_; - /// Identify those grid cells that are further subdivided by - /// an LGR. - /// - /// Writes to \c is_divided_. - /// - /// \param - void identifySubdividedCells(const ecl_grid_type* G); + /// Retrieve number of active cells in grid. + std::size_t numActiveCells() const; /// Compute linear index of global cell from explicit /// (I,J,K) tuple. @@ -374,6 +437,11 @@ namespace { /// Source cells for each Cartesian connection. OutCell outCell_; + /// Transmissibility field for purpose of on-demand flux + /// calculation if fluxes are not already available in dynamic + /// result set. + std::vector trans_; + /// Predicate for whether or not a particular result vector is /// defined on the grid's cells. /// @@ -415,6 +483,7 @@ namespace { void connectionData(const ::Opm::ECLResultData& src, const CartesianCells::Direction d, const std::string& vector, + const double unit, std::vector& x) const; /// Form complete name of directional result set vector from @@ -430,7 +499,7 @@ namespace { const CartesianCells::Direction d) const; /// Derive neighbourship relations on active cells in particular - /// Cartesian directions. + /// Cartesian directions and capture transmissibilty field. /// /// Writes to \c neigh_ and \c outCell_. /// @@ -487,6 +556,13 @@ ECL::getPVolVector(const ecl_grid_type* G, assert ((pvol.size() == nglob) && "Pore-volume must be provided for all global cells"); + + const auto pvol_unit = + getUnitSystem(init, gridID)->reservoirVolume(); + + for (auto& pv : pvol) { + pv = ::Opm::unit::convert::from(pv, pvol_unit); + } } return pvol; @@ -524,6 +600,18 @@ 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) + -> 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); + + return ::Opm::ECLUnits::createUnitSystem(ih[INTEHEAD_UNIT_INDEX]); +} + std::vector ECL::loadNNC(const ecl_grid_type* G, const ::Opm::ECLResultData& init) @@ -613,7 +701,7 @@ std::vector ECL::CartesianGridData::CartesianCells::activeGlobal() const { auto active = std::vector{}; - active.reserve(this->rsMap_.subset.size()); + active.reserve(this->numActiveCells()); for (const auto& id : this->rsMap_.subset) { active.push_back(id.glob); @@ -656,6 +744,48 @@ CartesianCells::scatterToGlobal(const std::vector& x) const return y; } +namespace { namespace ECL { + template + std::vector + CartesianGridData::CartesianCells:: + gatherToActive(const std::vector& x) const + { + const auto num_explicit_active = + static_cast(this->rsMap_.num_active); + + if (x.size() == num_explicit_active) { + // Input defined on explictly active cells (ACTNUM != 0). + // Extract subset of these. + + auto ax = std::vector{}; + ax.reserve(this->numActiveCells()); + + for (const auto& i : this->rsMap_.subset) { + ax.push_back(x[i.act]); + } + + return ax; + } + + if (x.size() == this->numGlobalCells()) { + // Input defined on all global cells. Extract active subset. + + auto ax = std::vector{}; + ax.reserve(this->numActiveCells()); + + for (const auto& i : this->rsMap_.subset) { + ax.push_back(x[i.glob]); + } + + return ax; + } + + // Input defined on neither explicitly active nor global cells. + // Possibly on all grid's NNCs. Let caller deal with this. + return x; + } +}} // namespace Anonymous::ECL + std::size_t ECL::CartesianGridData::CartesianCells::numGlobalCells() const { @@ -718,6 +848,12 @@ ECL::CartesianGridData::CartesianCells::isSubdivided(const int cellID) const return this->is_divided_[ix]; } +std::size_t +ECL::CartesianGridData::CartesianCells::numActiveCells() const +{ + return this->rsMap_.subset.size(); +} + std::size_t ECL::CartesianGridData:: CartesianCells::globIdx(const IndexTuple& ijk) const @@ -772,6 +908,7 @@ CartesianGridData(const ecl_grid_type* G, // Too large, but this is a quick estimate. this->neigh_.reserve(3 * (2 * this->numCells())); + this->trans_.reserve(3 * (1 * this->numCells())); for (const auto d : { CartesianCells::Direction::I , CartesianCells::Direction::J , @@ -781,6 +918,11 @@ CartesianGridData(const ecl_grid_type* G, } } +int ECL::CartesianGridData::gridID() const +{ + return this->gridID_; +} + std::size_t ECL::CartesianGridData::numCells() const { @@ -805,6 +947,12 @@ ECL::CartesianGridData::activePoreVolume() const return this->cells_.activePoreVolume(); } +const std::vector& +ECL::CartesianGridData::transmissibility() const +{ + return this->trans_; +} + int ECL::CartesianGridData::activeCell(const std::size_t globalCell) const { @@ -839,6 +987,22 @@ cellData(const ::Opm::ECLResultData& src, return this->cells_.scatterToGlobal(x); } +namespace { namespace ECL { + template + std::vector + CartesianGridData::activeCellData(const ::Opm::ECLResultData& src, + const std::string& vector) const + { + if (! this->haveCellData(src, vector)) { + return {}; + } + + auto x = src.keywordData(vector, this->gridID_); + + return this->cells_.gatherToActive(std::move(x)); + } +}} // namespace Anonymous::ECL + bool ECL::CartesianGridData:: haveCellData(const ::Opm::ECLResultData& src, @@ -871,7 +1035,8 @@ haveConnData(const ::Opm::ECLResultData& src, std::vector ECL::CartesianGridData:: connectionData(const ::Opm::ECLResultData& src, - const std::string& vector) const + const std::string& vector, + const double unit) const { if (! this->haveConnData(src, vector)) { return {}; @@ -883,7 +1048,7 @@ connectionData(const ::Opm::ECLResultData& src, CartesianCells::Direction::J , CartesianCells::Direction::K }) { - this->connectionData(src, d, vector, x); + this->connectionData(src, d, this->vectorName(vector, d), unit, x); } return x; @@ -894,10 +1059,10 @@ ECL::CartesianGridData:: connectionData(const ::Opm::ECLResultData& src, const CartesianCells::Direction d, const std::string& vector, + const double unit, std::vector& x) const { - const auto vname = this->vectorName(vector, d); - const auto v = this->cellData(src, vname); + const auto v = this->cellData(src, vector); const auto& cells = this->outCell_.find(d); @@ -905,7 +1070,7 @@ connectionData(const ::Opm::ECLResultData& src, "Direction must be I, J, or K"); for (const auto& cell : cells->second) { - x.push_back(v[cell]); + x.push_back(::Opm::unit::convert::from(v[cell], unit)); } } @@ -951,6 +1116,14 @@ deriveNeighbours(const std::vector& gcells, ? this->cellData(init, tran) : std::vector(this->cells_.numGlobalCells(), 1.0); + const auto trans_unit = + ECL::getUnitSystem(init, this->gridID_)->transmissibility(); + + auto SI_trans = [trans_unit](const double trans) + { + return ::Opm::unit::convert::from(trans, trans_unit); + }; + auto& ocell = this->outCell_[d]; ocell.reserve(gcells.size()); @@ -976,6 +1149,7 @@ deriveNeighbours(const std::vector& gcells, this->neigh_.push_back(other); ocell.push_back(globID); + this->trans_.push_back(SI_trans(T[globID])); } } } @@ -1032,28 +1206,57 @@ public: /// Retrieve number of active cells in graph. std::size_t numCells() const; - /// Retrive number of connections in graph. + /// Retrieve number of connections in graph. std::size_t numConnections() const; - /// Retrive neighbourship relations between active cells. + /// Retrieve the simulation scenario's active phases. + /// + /// Mostly useful to determine the set of \c PhaseIndex values for which + /// flux() may return non-zero values. + const std::vector& activePhases() const; + + /// Retrieve neighbourship relations between active cells. /// /// The \c i-th connection is between active cells \code /// neighbours()[2*i + 0] \endcode and \code neighbours()[2*i + 1] /// \endcode. std::vector neighbours() const; - /// Retrive static pore-volume values on active cells only. + /// Retrieve static pore-volume values on active cells only. /// /// Corresponds to the \c PORV vector in the INIT file, possibly /// restricted to those active cells for which the pore-volume is /// strictly positive. std::vector activePoreVolume() const; - const ::Opm::ECLResultData& rawResultData() const; + /// Retrieve static (background) transmissibility values on all + /// connections defined by \code neighbours() \endcode. + /// + /// Specifically, \code transmissibility()[i] \endcode is the + /// transmissibility of the connection between cells \code + /// neighbours()[2*i + 0] \endcode and \code neighbours()[2*i + 1] + /// \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; - /// Retrive phase flux on all connections defined by \code neighbours() + /// 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. /// /// \param[in] phase Canonical phase for which to retrive flux. @@ -1068,6 +1271,14 @@ public: std::vector flux(const PhaseIndex phase) const; + template + std::vector + rawLinearisedCellData(const std::string& vector) const; + + std::vector + linearisedCellData(const std::string& vector, + UnitConvention unit) const; + private: /// Collection of non-Cartesian neighbourship relations attributed to a /// particular ECL keyword set (i.e., one of NNC{1,2}, NNC{G,L}, NNCLL). @@ -1186,9 +1397,14 @@ private: /// \param[in] offset Start index into global linear number for all /// active grids. /// + /// \param[in] trans_unit Unit of measurement of transmissibility + /// field stored in result set. Used to convert values to the + /// strict SI conventions (i.e., rm^3). + /// /// \param[in] nnc Non-neighbouring connection from result set. void add(const std::vector& grids, const std::vector& offset, + const double trans_unit, const ecl_nnc_type& nnc); std::vector allCategories() const; @@ -1199,6 +1415,10 @@ private: /// Access all active non-neighbouring connections. const std::vector& getNeighbours() const; + /// Access transmissibility field of all active non-neighbouring + /// connections. Numerical values in strict SI units (rm^3). + const std::vector& transmissibility() const; + /// Retrieve all non-neighbouring connections of a particular /// category (i.e., pertaining to a particular set of keywords). /// @@ -1214,6 +1434,13 @@ private: /// in linear numbering of all model's active cells. std::vector neigh_; + /// Transmissibility of non-Cartesian (non-neighbouring) connections. + /// + /// Note that \code trans_[i] \endcode is the transmissibility of + /// the connection between cells \code neigh_[2*i + 0] \endcode and + /// \code neigh_[2*i + 1] \endcode. + std::vector trans_; + /// Collection of KeywordIndexMap keywords_; @@ -1221,7 +1448,9 @@ private: /// /// Simplifies implementation of ctor. /// - /// \param[in] cat Class + /// \param[in] cat Requested category of flux relation. + /// + /// \return Flux relation of type \p cat. FluxRelation makeRelation(const Category cat) const; /// Identify connection category from connection's grids. @@ -1284,6 +1513,11 @@ private: /// range \code [0 .. numCells()) \endcode). std::vector activeOffset_; + /// Set of active phases in result set. Derived from .INIT on the + /// assumption that the set of active phases does not change throughout + /// the simulation run. + std::vector activePhases_; + /// Current result set. std::unique_ptr src_; @@ -1299,6 +1533,11 @@ private: void defineNNCs(const ecl_grid_type* G, const ::Opm::ECLResultData& init); + /// Extract scenario's set of active phases. + /// + /// Writes to activePhases_. + void defineActivePhases(const ::Opm::ECLResultData& init); + /// Compute ECL vector basename for particular phase flux. /// /// \param[in] phase Canonical phase for which to derive ECL vector @@ -1312,15 +1551,29 @@ private: /// Extract flux values corresponding to particular result set vector /// for all identified non-neighbouring connections. /// + /// \tparam[in] GetFluxUnit Type of function object for computing the + /// grid-dependent flux unit. + /// /// \param[in] vector Result set vector prefix. Typically computed by /// method flowVector(). /// + /// \param[in] fluxUnit Function object for computing grid-dependent + /// flux units. Must support the syntax + /// \code + /// unit = fluxUnit(gridID) + /// \endcode + /// with 'gridID' being a non-negative \c int that identifies a + /// particular model grid (zero for the main grid and positive for + /// LGRs) and 'unit' a positive floating-point value. + /// /// \param[in,out] flux Numerical values of result set vector. On /// input, contains all values corresponding to all fully Cartesian /// connections across all active grids. On output additionally /// 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; }; @@ -1371,6 +1624,7 @@ void Opm::ECLGraph::Impl:: NNC::add(const std::vector& grid, const std::vector& offset, + const double trans_unit, const ecl_nnc_type& nnc) { if (! this->isViable(grid, nnc)) { @@ -1395,6 +1649,10 @@ NNC::add(const std::vector& grid, this->neigh_.push_back(o + c); } + // Capture transmissibility field to support on-demand flux calculations + // if flux fields are not output to the on-disk result set. + this->trans_.push_back(unit::convert::from(nnc.trans, trans_unit)); + const auto cat = this->classifyConnection(nnc.grid_nr1, nnc.grid_nr2); auto entry = NonNeighKeywordIndexSet::Map { @@ -1411,9 +1669,10 @@ NNC::add(const std::vector& grid, std::size_t Opm::ECLGraph::Impl::NNC::numConnections() const { - assert (this->neigh_.size() % 2 == 0); + assert ((this->neigh_.size() % 2) == 0); + assert ((this->neigh_.size() / 2) == this->trans_.size()); - return this->neigh_.size() / 2; + return this->trans_.size(); } const std::vector& @@ -1422,6 +1681,12 @@ Opm::ECLGraph::Impl::NNC::getNeighbours() const return this->neigh_; } +const std::vector& +Opm::ECLGraph::Impl::NNC::transmissibility() const +{ + return this->trans_; +} + const Opm::ECLGraph::Impl::NNC::FluxRelation& Opm::ECLGraph::Impl::NNC::getRelations(const Category& cat) const { @@ -1521,6 +1786,7 @@ Opm::ECLGraph::Impl::Impl(const Path& grid, const Path& init) } this->defineNNCs(G.get(), I); + this->defineActivePhases(I); } void @@ -1579,6 +1845,12 @@ Opm::ECLGraph::Impl::numConnections() const return nconn + this->nnc_.numConnections(); } +const std::vector& +Opm::ECLGraph::Impl::activePhases() const +{ + return this->activePhases_; +} + std::vector Opm::ECLGraph::Impl::neighbours() const { @@ -1625,6 +1897,34 @@ Opm::ECLGraph::Impl::activePoreVolume() const return pvol; } +std::vector +Opm::ECLGraph::Impl::transmissibility() const +{ + auto trans = std::vector{}; + + // Recall: this->numConnections() includes NNCs. + const auto totconn = this->numConnections(); + trans.reserve(totconn); + + for (const auto& G : this->grid_) { + const auto& Ti = G.transmissibility(); + + trans.insert(trans.end(), Ti.begin(), Ti.end()); + } + + if (this->nnc_.numConnections() > 0) { + const auto& tranNNC = this->nnc_.transmissibility(); + + trans.insert(trans.end(), tranNNC.begin(), tranNNC.end()); + } + + if (trans.size() < totconn) { + return {}; + } + + return trans; +} + const ::Opm::ECLResultData& Opm::ECLGraph::Impl::rawResultData() const { @@ -1640,6 +1940,11 @@ std::vector Opm::ECLGraph::Impl:: flux(const PhaseIndex phase) const { + auto fluxUnit = [this](const int gridID) + { + return ::ECL::getUnitSystem(*this->src_, gridID)->reservoirRate(); + }; + const auto vector = this->flowVector(phase); auto v = std::vector{}; @@ -1650,7 +1955,8 @@ flux(const PhaseIndex phase) const v.reserve(totconn); for (const auto& G : this->grid_) { - const auto& q = G.connectionData(*this->src_, vector); + const auto& q = + G.connectionData(*this->src_, vector, fluxUnit(G.gridID())); if (q.empty()) { // Flux vector invalid unless all grids provide this result @@ -1665,7 +1971,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, v); + this->fluxNNC(vector, std::move(fluxUnit), v); } if (v.size() < totconn) { @@ -1677,17 +1983,81 @@ flux(const PhaseIndex phase) const return v; } +namespace Opm { + + template + std::vector + ECLGraph::Impl::rawLinearisedCellData(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); + + x.insert(x.end(), std::begin(xi), std::end(xi)); + } + + if (x.size() != this->numCells()) { + return {}; + } + + return x; + } +} // namespace Opm + +std::vector +Opm::ECLGraph::Impl::linearisedCellData(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); + + if (xi.empty()) { continue; } + + // Note: Compensate for incrementing Grid ID above. + const auto usys = + ECL::getUnitSystem(*this->src_, G.gridID()); + + // Note: 'usys' (generally, std::unique_ptr<>) does not support + // regular PMF syntax (i.e., operator->*()). + const auto vector_unit = ((*usys).*unit)(); + + std::transform(std::begin(xi), std::end(xi), + std::back_inserter(x), + [vector_unit](const double value) + { + return ::Opm::unit::convert::from(value, vector_unit); + }); + } + + if (x.size() != this->numCells()) { + return {}; + } + + return x; +} + void Opm::ECLGraph::Impl::defineNNCs(const ecl_grid_type* G, const ::Opm::ECLResultData& init) { + // Assume all transmissibilites in the result set follow the same unit + // conventions. + + const auto trans_unit = + ECL::getUnitSystem(init, 0)->transmissibility(); + for (const auto& nnc : ECL::loadNNC(G, init)) { - this->nnc_.add(this->grid_, this->activeOffset_, nnc); + 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 { auto v = std::vector(this->nnc_.numConnections(), 0.0); @@ -1697,13 +2067,11 @@ Opm::ECLGraph::Impl::fluxNNC(const std::string& vector, const auto& rel = this->nnc_.getRelations(cat); const auto fluxID = rel.makeKeyword(vector); - auto gridID = 0; for (const auto& G : this->grid_) { - const auto& iset = rel.indexSet().getGridCollection(gridID); + const auto gridID = G.gridID(); - // Must increment grid ID irrespective of early break of - // iteration. - gridID += 1; + const auto& iset = + rel.indexSet().getGridCollection(gridID); if (iset.empty()) { // No NNCs for this category in this grid. Skip. @@ -1719,13 +2087,16 @@ Opm::ECLGraph::Impl::fluxNNC(const std::string& vector, continue; } + const auto flux_unit = fluxUnit(gridID); + // Data fully available for (category,gridID). Assign // approriate subset of NNC flux vector. for (const auto& ix : iset) { assert (ix.neighIdx < v.size()); assert (ix.kwIdx < q.size()); - v[ix.neighIdx] = q[ix.kwIdx]; + v[ix.neighIdx] = + unit::convert::from(q[ix.kwIdx], flux_unit); assigned[ix.neighIdx] = true; } @@ -1745,6 +2116,32 @@ Opm::ECLGraph::Impl::fluxNNC(const std::string& vector, } } +void +Opm::ECLGraph::Impl:: +defineActivePhases(const ::Opm::ECLResultData& init) +{ + const auto ih = + init.keywordData(INTEHEAD_KW, ECL_GRID_MAINGRID_LGR_NR); + + const auto phaseMask = + static_cast(ih[INTEHEAD_PHASE_INDEX]); + + using Check = std::pair; + + const auto allPhases = std::vector { + { PhaseIndex::Aqua , (1u << 1) }, + { PhaseIndex::Liquid, (1u << 0) }, + { PhaseIndex::Vapour, (1u << 2) }, + }; + + this->activePhases_.clear(); + for (const auto& phase : allPhases) { + if ((phase.second & phaseMask) != 0) { + this->activePhases_.push_back(phase.first); + } + } +} + std::string Opm::ECLGraph::Impl:: flowVector(const PhaseIndex phase) const @@ -1834,6 +2231,12 @@ Opm::ECLGraph::numConnections() const return this->pImpl_->numConnections(); } +const std::vector& +Opm::ECLGraph::activePhases() const +{ + return this->pImpl_->activePhases(); +} + std::vector Opm::ECLGraph::neighbours() const { return this->pImpl_->neighbours(); @@ -1844,6 +2247,11 @@ std::vector Opm::ECLGraph::poreVolume() const return this->pImpl_->activePoreVolume(); } +std::vector Opm::ECLGraph::transmissibility() const +{ + return this->pImpl_->transmissibility(); +} + bool Opm::ECLGraph::selectReportStep(const int rptstep) const { return this->pImpl_->selectReportStep(rptstep); @@ -1861,3 +2269,28 @@ flux(const PhaseIndex phase) const { return this->pImpl_->flux(phase); } + +namespace Opm { + + template + std::vector + ECLGraph::rawLinearisedCellData(const std::string& vector) const + { + return this->pImpl_->rawLinearisedCellData(vector); + } + + // Explicit instantiations for the element types we care about. + template std::vector + ECLGraph::rawLinearisedCellData(const std::string& vector) const; + + template std::vector + ECLGraph::rawLinearisedCellData(const std::string& vector) const; + +} // namespace Opm + +std::vector +Opm::ECLGraph::linearisedCellData(const std::string& vector, + UnitConvention unit) const +{ + return this->pImpl_->linearisedCellData(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 c2fc6462a9..62204e5722 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 @@ -22,6 +22,7 @@ #define OPM_ECLGRAPH_HEADER_INCLUDED #include +#include #include #include @@ -46,6 +47,9 @@ namespace Opm { public: using Path = boost::filesystem::path; + /// Enum for indicating requested phase from the flux() method. + enum class PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 }; + /// Disabled default constructor. ECLGraph() = delete; @@ -124,6 +128,12 @@ namespace Opm { /// Retrieve number of connections in graph. std::size_t numConnections() const; + /// Retrieve the simulation scenario's active phases. + /// + /// Mostly useful to determine the set of \c PhaseIndex values for + /// which flux() will return non-zero values if data available. + const std::vector& activePhases() const; + /// Retrieve neighbourship relations between active cells. /// /// The \c i-th connection is between active cells \code @@ -135,9 +145,18 @@ namespace Opm { /// /// Corresponds to the \c PORV vector in the INIT file, possibly /// restricted to those active cells for which the pore-volume is - /// strictly positive. + /// strictly positive. Numerical values in SI units (rm^3). std::vector poreVolume() const; + /// Retrieve static (background) transmissibility values on all + /// connections defined by \code neighbours() \endcode. + /// + /// Specifically, \code transmissibility()[i] \endcode is the + /// transmissibility of the connection between cells \code + /// neighbours()[2*i + 0] \endcode and \code neighbours()[2*i + 1] + /// \endcode. + std::vector transmissibility() const; + /// Restrict dynamic result set data to single report step. /// /// This method must be called before calling either flux() or @@ -156,21 +175,57 @@ namespace Opm { /// to method selectReportStep(). const ::Opm::ECLResultData& rawResultData() const; - /// Enum for indicating requested phase from the flux() method. - enum class PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 }; - - /// Retrive phase flux on all connections defined by \code + /// Retrieve phase flux on all connections defined by \code /// neighbours() \endcode. /// - /// \param[in] phase Canonical phase for which to retrive flux. + /// \param[in] phase Canonical phase for which to retrieve flux. /// /// \return Flux values corresponding to selected phase. Empty if /// unavailable in the result set (e.g., when querying the gas /// flux in an oil/water system or if no flux values at all were - /// output to the restart file). + /// output to the restart file). Numerical values in SI units + /// (rm^3/s). std::vector flux(const PhaseIndex phase) const; + /// Retrieve result set vector from current view (e.g., particular + /// report step) linearised on active cells. + /// + /// \tparam T Element type of result set vector. + /// + /// \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; + + /// Convenience type alias for \c UnitSystem PMFs (pointer to member + /// function). + typedef double (ECLUnits::UnitSystem::*UnitConvention)() 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("PRESSURE", &ECLUnits::UnitSystem::pressure); + /// \endcode + /// + /// \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; + private: /// Implementation class. class Impl; 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 new file mode 100644 index 0000000000..4a3b77b1c1 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.cpp @@ -0,0 +1,209 @@ +/* + Copyright 2017 SINTEF ICT, Applied Mathematics. + Copyright 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 + +#include + +#include +#include +#include + +#include + +namespace Opm { namespace ECLUnits { + namespace Impl { + ert_ecl_unit_enum getUnitConvention(const int usys); + + template + class USys; + + template <> + class USys : public ::Opm::ECLUnits::UnitSystem + { + public: + virtual double pressure() const override + { + return Metric::Pressure; + } + + virtual double reservoirRate() const override + { + return Metric::ReservoirVolume / Metric::Time; + } + + virtual double reservoirVolume() const override + { + return Metric::ReservoirVolume; + } + + virtual double time() const override + { + return Metric::Time; + } + + virtual double transmissibility() const override + { + return Metric::Transmissibility; + } + }; + + template <> + class USys : public ::Opm::ECLUnits::UnitSystem + { + public: + virtual double pressure() const override + { + return Field::Pressure; + } + + virtual double reservoirRate() const override + { + return Field::ReservoirVolume / Field::Time; + } + + virtual double reservoirVolume() const override + { + return Field::ReservoirVolume; + } + + virtual double time() const override + { + return Field::Time; + } + + virtual double transmissibility() const override + { + return Field::Transmissibility; + } + }; + + template <> + class USys : public ::Opm::ECLUnits::UnitSystem + { + public: + virtual double pressure() const override + { + return Lab::Pressure; + } + + virtual double reservoirRate() const override + { + return Lab::ReservoirVolume / Lab::Time; + } + + virtual double reservoirVolume() const override + { + return Lab::ReservoirVolume; + } + + virtual double time() const override + { + return Lab::Time; + } + + virtual double transmissibility() const override + { + return Lab::Transmissibility; + } + }; + + template <> + class USys : public ::Opm::ECLUnits::UnitSystem + { + public: + virtual double pressure() const override + { + return unit::atm; + } + + virtual double reservoirRate() const override + { + using namespace prefix; + using namespace unit; + + return cubic(meter) / day; + } + + virtual double reservoirVolume() const override + { + using namespace prefix; + using namespace unit; + + return cubic(meter); + } + + virtual double time() const override + { + return unit::day; + } + + virtual double transmissibility() const override + { + using namespace prefix; + using namespace unit; + + return centi*Poise * cubic(meter) / (day * atm); + } + }; + } // namespace Impl +}} // namespace Opm::ECLUnits + +ert_ecl_unit_enum +Opm::ECLUnits::Impl::getUnitConvention(const int usys) +{ + switch (usys) { + case 1: return ECL_METRIC_UNITS; + case 2: return ECL_FIELD_UNITS; + case 3: return ECL_LAB_UNITS; + case 4: return ECL_PVT_M_UNITS; + } + + throw std::runtime_error("Unsupported Unit Convention: " + + std::to_string(usys)); +} + +std::unique_ptr +Opm::ECLUnits::createUnitSystem(const int usys) +{ + using UPtr = + std::unique_ptr; + + switch (Impl::getUnitConvention(usys)) { + case ECL_METRIC_UNITS: + return UPtr{ new Impl::USys{} }; + + case ECL_FIELD_UNITS: + return UPtr{ new Impl::USys{} }; + + case ECL_LAB_UNITS: + return UPtr{ new Impl::USys{} }; + + case ECL_PVT_M_UNITS: + return UPtr{ new Impl::USys{} }; + } + + throw std::runtime_error("Unsupported Unit Convention: " + + std::to_string(usys)); +} 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 new file mode 100644 index 0000000000..ee11974fe0 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLUnitHandling.hpp @@ -0,0 +1,45 @@ +/* + Copyright 2017 SINTEF ICT, Applied Mathematics. + Copyright 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 . +*/ + +#ifndef OPM_ECLUNITHANDLING_HEADER_INCLUDED +#define OPM_ECLUNITHANDLING_HEADER_INCLUDED + +#include + +namespace Opm { + + namespace ECLUnits { + + struct UnitSystem + { + virtual double pressure() const = 0; + virtual double reservoirRate() const = 0; + virtual double reservoirVolume() const = 0; + virtual double time() const = 0; + virtual double transmissibility() const = 0; + }; + + std::unique_ptr + createUnitSystem(const int usys); + + } // namespace ECLUnits +} // namespace Opm + +#endif // OPM_ECLUNITHANDLING_HEADER_INCLUDED 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 630f17ca89..20b611a4a0 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 @@ -24,6 +24,7 @@ #include #include +#include #include #include #include @@ -108,15 +109,7 @@ namespace Opm // Reservoir rate units from code used in INTEHEAD. double resRateUnit(const int unit_code) { - using prefix::centi; - using namespace unit; - - switch (unit_code) { - case INTEHEAD::Metric: return cubic(meter)/day; // m^3/day - case INTEHEAD::Field: return stb/day; // stb/day - case INTEHEAD::Lab: return cubic(centi*meter)/hour; // (cm)^3/h - default: throw std::runtime_error("Unknown unit code from INTEHEAD: " + std::to_string(unit_code)); - } + return ECLUnits::createUnitSystem(unit_code)->reservoirRate(); } 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 new file mode 100644 index 0000000000..2d9051cf2e --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runAcceptanceTest.cpp @@ -0,0 +1,550 @@ +/* + Copyright 2017 SINTEF ICT, Applied Mathematics. + Copyright 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 . +*/ + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace StringUtils { + namespace { + std::string trim(const std::string& s) + { + const auto anchor_ws = + boost::regex(R"~~(^\s+([^\s]+)\s+$)~~"); + + auto m = boost::smatch{}; + + if (boost::regex_match(s, m, anchor_ws)) { + return m[1]; + } + + return s; + } + + std::vector split(const std::string& s) + { + if (s.empty()) { + // Single element vector whose only element is the empty + // string. + return { "" }; + } + + const auto sep = boost::regex(R"~~([\s,;.|]+)~~"); + + using TI = boost::sregex_token_iterator; + + // vector(begin, end) + // + // Range is every substring (i.e., token) in input string 's' + // that does NOT match 'sep'. + return{ TI(s.begin(), s.end(), sep, -1), TI{} }; + } + } // namespace Anonymous + + template + struct StringTo; + + template <> + struct StringTo + { + static int value(const std::string& s); + }; + + template <> + struct StringTo + { + static double value(const std::string& s); + }; + + template <> + struct StringTo + { + static std::string value(const std::string& s); + }; + + int StringTo::value(const std::string& s) + { + return std::stoi(s); + } + + double StringTo::value(const std::string& s) + { + return std::stod(s); + } + + std::string StringTo::value(const std::string& s) + { + return trim(s); + } + + namespace VectorValue { + template + std::vector get(const std::string& s, std::true_type) + { + return split(s); + } + + template + std::vector get(const std::string& s, std::false_type) + { + const auto tokens = split(s); + + auto ret = std::vector{}; + ret.reserve(tokens.size()); + for (const auto& token : tokens) { + ret.push_back(StringTo::value(token)); + } + + return ret; + } + + template + std::vector get(const std::string& s) + { + return get(s, typename std::is_same::type()); + } + } // namespace VectorValue +} // namespace StringUtils + +namespace { + struct PoreVolume + { + std::vector data; + }; + + class VectorDifference + { + public: + using Vector = std::vector; + using size_type = Vector::size_type; + + VectorDifference(const Vector& x, const Vector& y) + : x_(x), y_(y) + { + if (x_.size() != y_.size()) { + std::ostringstream os; + + os << "Incompatible Array Sizes: Expected 2x" + << x_.size() << ", but got (" + << x_.size() << ", " << y_.size() << ')'; + + throw std::domain_error(os.str()); + } + } + + size_type size() const + { + return x_.size(); + } + + bool empty() const + { + return this->size() == 0; + } + + double operator[](const size_type i) const + { + return x_[i] - y_[i]; + } + + private: + const Vector& x_; + const Vector& y_; + }; + + template + class VectorRatio + { + public: + using size_type = typename std::decay< + decltype(std::declval()[0]) + >::type; + + VectorRatio(const Vector1& x, const Vector2& y) + : x_(x), y_(y) + { + if (x_.size() != y.size()) { + std::ostringstream os; + + os << "Incompatible Array Sizes: Expected 2x" + << x_.size() << ", but got (" + << x_.size() << ", " << y_.size() << ')'; + + throw std::domain_error(os.str()); + } + } + + size_type size() const + { + return x_.size(); + } + + bool empty() const + { + return x_.empty(); + } + + double operator[](const size_type i) const + { + return x_[i] / y_[i]; + } + + private: + const Vector1& x_; + const Vector2& y_; + }; + + struct ErrorMeasurement + { + double volume; + double inf; + }; + + struct ErrorTolerance + { + double absolute; + double relative; + }; + + struct AggregateErrors + { + std::vector absolute; + std::vector relative; + }; + + struct ReferenceToF + { + std::vector forward; + std::vector reverse; + }; + + template + double volumeMetric(const PoreVolume& pv, + const FieldVariable& x) + { + if (x.size() != pv.data.size()) { + std::ostringstream os; + + os << "Incompatible Array Sizes: Expected 2x" + << pv.data.size() << ", but got (" + << pv.data.size() << ", " << x.size() << ')'; + + throw std::domain_error(os.str()); + } + + auto num = 0.0; + auto den = 0.0; + + for (decltype(pv.data.size()) + i = 0, n = pv.data.size(); i < n; ++i) + { + num += std::abs(x[i]) * pv.data[i]; + den += pv.data[i]; + } + + return num / den; + } + + template + double pointMetric(const FieldVariable& x) + { + static_assert(std::is_same::type, double>::value, + "Field Variable Value Type Must be 'double'"); + + if (x.empty()) { + return 0; + } + + auto max = 0*x[0] - 1; + + for (decltype(x.size()) + i = 0, n = x.size(); i < n; ++i) + { + const auto t = std::abs(x[i]); + + if (t > max) { + max = t; + } + } + + return max; + } + + std::vector + availableReportSteps(const example::FilePaths& paths) + { + using FilePtr = ::ERT:: + ert_unique_ptr; + + const auto rsspec_fn = example:: + deriveFileName(paths.grid, { ".RSSPEC", ".FRSSPEC" }); + + // Read-only, keep open between requests + const auto open_flags = 0; + + auto rsspec = FilePtr{ + ecl_file_open(rsspec_fn.generic_string().c_str(), open_flags) + }; + + auto* globView = ecl_file_get_global_view(rsspec.get()); + + const auto* ITIME_kw = "ITIME"; + const auto n = ecl_file_view_get_num_named_kw(globView, ITIME_kw); + + auto steps = std::vector(n); + + for (auto i = 0*n; i < n; ++i) { + const auto* itime = + ecl_file_view_iget_named_kw(globView, ITIME_kw, i); + + const auto* itime_data = + static_cast(ecl_kw_iget_ptr(itime, 0)); + + steps[i] = itime_data[0]; + } + + return steps; + } + + ErrorTolerance + testTolerances(const ::Opm::parameter::ParameterGroup& param) + { + const auto atol = param.getDefault("atol", 1.0e-8); + const auto rtol = param.getDefault("rtol", 5.0e-12); + + return ErrorTolerance{ atol, rtol }; + } + + int numDigits(const std::vector& steps) + { + if (steps.empty()) { + return 1; + } + + const auto m = + *std::max_element(std::begin(steps), std::end(steps)); + + if (m == 0) { + return 1; + } + + assert (m > 0); + + return std::floor(std::log10(static_cast(m))) + 1; + } + + ReferenceToF + loadReference(const ::Opm::parameter::ParameterGroup& param, + const int step, + const int nDigits) + { + namespace fs = boost::filesystem; + + using VRef = std::reference_wrapper>; + + auto fname = fs::path(param.get("ref-dir")); + { + std::ostringstream os; + + os << "tof-" << std::setw(nDigits) << std::setfill('0') + << step << ".txt"; + + fname /= os.str(); + } + + fs::ifstream input(fname); + + if (! input) { + std::ostringstream os; + + os << "Unable to Open Reference Data File " + << fname.filename(); + + throw std::domain_error(os.str()); + } + + auto tof = ReferenceToF{}; + + auto ref = std::array{{ std::ref(tof.forward) , + std::ref(tof.reverse) }}; + + { + auto i = static_cast(0); + auto t = 0.0; + + while (input >> t) { + ref[i].get().push_back(t); + + i = (i + 1) % 2; + } + } + + if (tof.forward.size() != tof.reverse.size()) { + std::ostringstream os; + + os << "Unable to Read Consistent ToF Reference Data From " + << fname.filename(); + + throw std::out_of_range(os.str()); + } + + return tof; + } + + void computeErrors(const PoreVolume& pv, + const std::vector& ref, + const ::Opm::FlowDiagnostics::Solution& fd, + AggregateErrors& E) + { + const auto tof = fd.timeOfFlight(); + const auto diff = VectorDifference(tof, ref); // tof - ref + + using Vector1 = std::decay::type; + using Vector2 = std::decay::type; + using Ratio = VectorRatio; + + const auto rat = Ratio(diff, ref); // (tof - ref) / ref + + auto abs = ErrorMeasurement{}; + { + abs.volume = volumeMetric(pv, diff); + abs.inf = pointMetric ( diff); + } + + auto rel = ErrorMeasurement{}; + { + rel.volume = volumeMetric(pv, rat); + rel.inf = pointMetric ( rat); + } + + E.absolute.push_back(std::move(abs)); + E.relative.push_back(std::move(rel)); + } + + std::array + sampleDifferences(example::Setup&& setup, + const std::vector& steps) + { + const auto start = + std::vector{}; + + const auto nDigits = numDigits(steps); + + const auto pv = PoreVolume{ setup.graph.poreVolume() }; + + auto E = std::array{}; + + for (const auto& step : steps) { + if (step == 0) { + // Ignore initial condition + continue; + } + + if (! setup.selectReportStep(step)) { + continue; + } + + const auto ref = loadReference(setup.param, step, nDigits); + + { + const auto fwd = setup.toolbox + .computeInjectionDiagnostics(start); + + computeErrors(pv, ref.forward, fwd.fd, E[0]); + } + + { + const auto rev = setup.toolbox + .computeProductionDiagnostics(start); + + computeErrors(pv, ref.reverse, rev.fd, E[1]); + } + } + + return E; + } + + bool errorAcceptable(const std::vector& E, + const double tol) + { + return std::accumulate(std::begin(E), std::end(E), true, + [tol](const bool ok, const ErrorMeasurement& e) + { + // Fine if at least one of .volume or .inf <= tol. + return ok && ! ((e.volume > tol) && (e.inf > tol)); + }); + } + + bool everythingFine(const AggregateErrors& E, + const ErrorTolerance& tol) + { + return errorAcceptable(E.absolute, tol.absolute) + && errorAcceptable(E.relative, tol.relative); + } +} // namespace Anonymous + +int main(int argc, char* argv[]) +try { + auto setup = example::Setup(argc, argv); + + const auto tol = testTolerances(setup.param); + const auto steps = availableReportSteps(setup.file_paths); + + const auto E = sampleDifferences(std::move(setup), steps); + const auto ok = + everythingFine(E[0], tol) && everythingFine(E[1], tol); + + std::cout << (ok ? "OK" : "FAIL") << '\n'; + + if (! ok) { + return EXIT_FAILURE; + } +} +catch (const std::exception& e) { + std::cerr << "Caught Exception: " << e.what() << '\n'; + + return EXIT_FAILURE; +} 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 new file mode 100644 index 0000000000..89500d9b8e --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runLinearisedCellDataTest.cpp @@ -0,0 +1,587 @@ +/* + Copyright 2017 SINTEF ICT, Applied Mathematics. + Copyright 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 . +*/ + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace StringUtils { + namespace { + std::string trim(const std::string& s) + { + const auto anchor_ws = + boost::regex(R"~~(^\s+([^\s]+)\s+$)~~"); + + auto m = boost::smatch{}; + + if (boost::regex_match(s, m, anchor_ws)) { + return m[1]; + } + + return s; + } + + std::vector split(const std::string& s) + { + if (s.empty()) { + // Single element vector whose only element is the empty + // string. + return { "" }; + } + + const auto sep = boost::regex(R"~~([\s,;.|]+)~~"); + + using TI = boost::sregex_token_iterator; + + // vector(begin, end) + // + // Range is every substring (i.e., token) in input string 's' + // that does NOT match 'sep'. + return{ TI(s.begin(), s.end(), sep, -1), TI{} }; + } + } // namespace Anonymous + + template + struct StringTo; + + template <> + struct StringTo + { + static int value(const std::string& s); + }; + + template <> + struct StringTo + { + static double value(const std::string& s); + }; + + template <> + struct StringTo + { + static std::string value(const std::string& s); + }; + + int StringTo::value(const std::string& s) + { + return std::stoi(s); + } + + double StringTo::value(const std::string& s) + { + return std::stod(s); + } + + std::string StringTo::value(const std::string& s) + { + return trim(s); + } + + namespace VectorValue { + template + std::vector get(const std::string& s, std::true_type) + { + return split(s); + } + + template + std::vector get(const std::string& s, std::false_type) + { + const auto tokens = split(s); + + auto ret = std::vector{}; + ret.reserve(tokens.size()); + for (const auto& token : tokens) { + ret.push_back(StringTo::value(token)); + } + + return ret; + } + + template + std::vector get(const std::string& s) + { + return get(s, typename std::is_same::type()); + } + } // namespace VectorValue +} // namespace StringUtils + +namespace { + struct Reference + { + std::vector data; + }; + + struct Calculated + { + std::vector data; + }; + + class VectorUnits + { + private: + using USys = ::Opm::ECLUnits::UnitSystem; + + public: + using UnitConvention = ::Opm::ECLGraph::UnitConvention; + + VectorUnits() + : units_({ { "pressure", &USys::pressure } }) + { + } + + UnitConvention getUnit(const std::string& vector) const + { + auto p = units_.find(vector); + + if (p == units_.end()) { + std::ostringstream os; + + os << "Unsupported Vector Quantity '" << vector << '\''; + + throw std::domain_error(os.str()); + } + + return p->second; + } + + private: + std::map units_; + }; + + class VectorDifference + { + public: + using Vector = std::vector; + using size_type = Vector::size_type; + + VectorDifference(const Vector& x, const Vector& y) + : x_(x), y_(y) + { + if (x_.size() != y_.size()) { + std::ostringstream os; + + os << "Incompatible Array Sizes: Expected 2x" + << x_.size() << ", but got (" + << x_.size() << ", " << y_.size() << ')'; + + throw std::domain_error(os.str()); + } + } + + size_type size() const + { + return x_.size(); + } + + bool empty() const + { + return this->size() == 0; + } + + double operator[](const size_type i) const + { + return x_[i] - y_[i]; + } + + private: + const Vector& x_; + const Vector& y_; + }; + + template + class VectorRatio + { + public: + using size_type = typename std::decay< + decltype(std::declval()[0]) + >::type; + + VectorRatio(const Vector1& x, const Vector2& y) + : x_(x), y_(y) + { + if (x_.size() != y.size()) { + std::ostringstream os; + + os << "Incompatible Array Sizes: Expected 2x" + << x_.size() << ", but got (" + << x_.size() << ", " << y_.size() << ')'; + + throw std::domain_error(os.str()); + } + } + + size_type size() const + { + return x_.size(); + } + + bool empty() const + { + return x_.empty(); + } + + double operator[](const size_type i) const + { + return x_[i] / y_[i]; + } + + private: + const Vector1& x_; + const Vector2& y_; + }; + + struct ErrorMeasurement + { + double volume; + double inf; + }; + + struct ErrorTolerance + { + double absolute; + double relative; + }; + + struct AggregateErrors + { + std::vector absolute; + std::vector relative; + }; + + struct ReferenceSolution + { + std::vector raw; + std::vector SI; + }; + + template + double volumeMetric(const FieldVariable& x) + { + auto result = 0.0; + + for (decltype(x.size()) + i = 0, n = x.size(); i < n; ++i) + { + const auto m = std::abs(x[i]); + result += m * m; + } + + return std::sqrt(result / x.size()); + } + + template + double pointMetric(const FieldVariable& x) + { + static_assert(std::is_same::type, double>::value, + "Field Variable Value Type Must be 'double'"); + + if (x.empty()) { + return 0; + } + + auto max = 0*x[0] - 1; + + for (decltype(x.size()) + i = 0, n = x.size(); i < n; ++i) + { + const auto t = std::abs(x[i]); + + if (t > max) { + max = t; + } + } + + return max; + } + + std::vector + availableReportSteps(const example::FilePaths& paths) + { + using FilePtr = ::ERT:: + ert_unique_ptr; + + const auto rsspec_fn = example:: + deriveFileName(paths.grid, { ".RSSPEC", ".FRSSPEC" }); + + // Read-only, keep open between requests + const auto open_flags = 0; + + auto rsspec = FilePtr{ + ecl_file_open(rsspec_fn.generic_string().c_str(), open_flags) + }; + + auto* globView = ecl_file_get_global_view(rsspec.get()); + + const auto* ITIME_kw = "ITIME"; + const auto n = ecl_file_view_get_num_named_kw(globView, ITIME_kw); + + auto steps = std::vector(n); + + for (auto i = 0*n; i < n; ++i) { + const auto* itime = + ecl_file_view_iget_named_kw(globView, ITIME_kw, i); + + const auto* itime_data = + static_cast(ecl_kw_iget_ptr(itime, 0)); + + steps[i] = itime_data[0]; + } + + return steps; + } + + ErrorTolerance + testTolerances(const ::Opm::parameter::ParameterGroup& param) + { + const auto atol = param.getDefault("atol", 1.0e-8); + const auto rtol = param.getDefault("rtol", 5.0e-12); + + return ErrorTolerance{ atol, rtol }; + } + + std::vector + testQuantities(const ::Opm::parameter::ParameterGroup& param) + { + return StringUtils::VectorValue:: + get(param.get("quant")); + } + + int numDigits(const std::vector& steps) + { + if (steps.empty()) { + return 1; + } + + const auto m = + *std::max_element(std::begin(steps), std::end(steps)); + + if (m == 0) { + return 1; + } + + assert (m > 0); + + return std::floor(std::log10(static_cast(m))) + 1; + } + + ReferenceSolution + loadReference(const ::Opm::parameter::ParameterGroup& param, + const std::string& quant, + const int step, + const int nDigits) + { + namespace fs = boost::filesystem; + + using VRef = std::reference_wrapper>; + + auto x = ReferenceSolution{}; + auto ref = std::array{{ std::ref(x.raw) , + std::ref(x.SI ) }}; + + auto i = 0; + + for (const auto* q : { "raw", "SI" }) { + auto fname = fs::path(param.get("ref-dir")) + / boost::algorithm::to_lower_copy(quant); + { + std::ostringstream os; + + os << q << '-' + << std::setw(nDigits) << std::setfill('0') + << step << ".txt"; + + fname /= os.str(); + } + + fs::ifstream input(fname); + + if (input) { + ref[i].get().assign(std::istream_iterator(input), + std::istream_iterator()); + } + + i += 1; + } + + if (x.raw.size() != x.SI.size()) { + std::ostringstream os; + + os << "Unable to Read Consistent Reference Data From '" + << param.get("ref-dir") << "' In Step " + << step; + + throw std::out_of_range(os.str()); + } + + return x; + } + + void computeErrors(const Reference& ref, + const Calculated& calc, + AggregateErrors& E) + { + const auto diff = + VectorDifference(calc.data, ref.data); // calc - ref + + using Vector1 = std::decay::type; + using Vector2 = std::decay::type; + using Ratio = VectorRatio; + + const auto rat = Ratio(diff, ref.data); // (tof - ref) / ref + + auto abs = ErrorMeasurement{}; + { + abs.volume = volumeMetric(diff); + abs.inf = pointMetric (diff); + } + + auto rel = ErrorMeasurement{}; + { + rel.volume = volumeMetric(rat); + rel.inf = pointMetric (rat); + } + + E.absolute.push_back(std::move(abs)); + E.relative.push_back(std::move(rel)); + } + + std::array + sampleDifferences(const ::Opm::ECLGraph& graph, + const ::Opm::parameter::ParameterGroup& param, + const std::string& quant, + const std::vector& steps) + { + const auto ECLquant = boost::algorithm::to_upper_copy(quant); + + auto unit = VectorUnits() + .getUnit(boost::algorithm::to_lower_copy(quant)); + + const auto start = + std::vector{}; + + const auto nDigits = numDigits(steps); + + auto E = std::array{}; + + for (const auto& step : steps) { + if (! graph.selectReportStep(step)) { + continue; + } + + const auto ref = loadReference(param, quant, step, nDigits); + + { + const auto raw = Calculated { + graph.rawLinearisedCellData(ECLquant) + }; + + computeErrors(Reference{ ref.raw }, raw, E[0]); + } + + { + const auto SI = Calculated { + graph.linearisedCellData(ECLquant, unit) + }; + + computeErrors(Reference{ ref.SI }, SI, E[1]); + } + } + + return E; + } + + bool errorAcceptable(const std::vector& E, + const double tol) + { + return std::accumulate(std::begin(E), std::end(E), true, + [tol](const bool ok, const ErrorMeasurement& e) + { + // Fine if at least one of .volume or .inf <= tol. + return ok && ! ((e.volume > tol) && (e.inf > tol)); + }); + } + + bool everythingFine(const AggregateErrors& E, + const ErrorTolerance& tol) + { + return errorAcceptable(E.absolute, tol.absolute) + && errorAcceptable(E.relative, tol.relative); + } +} // namespace Anonymous + +int main(int argc, char* argv[]) +try { + const auto prm = example::initParam(argc, argv); + const auto pth = example::FilePaths(prm); + const auto tol = testTolerances(prm); + + 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 ok = + everythingFine(E[0], tol) && everythingFine(E[1], tol); + + std::cout << quant << ": " << (ok ? "OK" : "FAIL") << '\n'; + + all_ok = all_ok && ok; + } + + if (! all_ok) { + return EXIT_FAILURE; + } +} +catch (const std::exception& e) { + std::cerr << "Caught Exception: " << e.what() << '\n'; + + return EXIT_FAILURE; +} 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 new file mode 100644 index 0000000000..e705cfede4 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runTransTest.cpp @@ -0,0 +1,229 @@ +/* + Copyright 2017 SINTEF ICT, Applied Mathematics. + Copyright 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 . +*/ + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace { + class VectorDifference + { + public: + using Vector = std::vector; + using size_type = Vector::size_type; + + VectorDifference(const Vector& x, const Vector& y) + : x_(x), y_(y) + { + if (x_.size() != y_.size()) { + std::ostringstream os; + + os << "Incompatible Array Sizes: Expected 2x" + << x_.size() << ", but got (" + << x_.size() << ", " << y_.size() << ')'; + + throw std::domain_error(os.str()); + } + } + + size_type size() const + { + return x_.size(); + } + + bool empty() const + { + return this->size() == 0; + } + + double operator[](const size_type i) const + { + return x_[i] - y_[i]; + } + + private: + const Vector& x_; + const Vector& y_; + }; + + template + class VectorRatio + { + public: + using size_type = typename std::decay< + decltype(std::declval()[0]) + >::type; + + VectorRatio(const Vector1& x, const Vector2& y) + : x_(x), y_(y) + { + if (x_.size() != y.size()) { + std::ostringstream os; + + os << "Incompatible Array Sizes: Expected 2x" + << x_.size() << ", but got (" + << x_.size() << ", " << y_.size() << ')'; + + throw std::domain_error(os.str()); + } + } + + size_type size() const + { + return x_.size(); + } + + bool empty() const + { + return x_.empty(); + } + + double operator[](const size_type i) const + { + return x_[i] / y_[i]; + } + + private: + const Vector1& x_; + const Vector2& y_; + }; + + struct ErrorTolerance + { + double absolute; + double relative; + }; + + template + double pointMetric(const FieldVariable& x) + { + static_assert(std::is_same::type, double>::value, + "Field Variable Value Type Must be 'double'"); + + if (x.empty()) { + return 0; + } + + auto max = 0*x[0] - 1; + + for (decltype(x.size()) + i = 0, n = x.size(); i < n; ++i) + { + const auto t = std::abs(x[i]); + + if (t > max) { + max = t; + } + } + + return max; + } + + ErrorTolerance + testTolerances(const ::Opm::parameter::ParameterGroup& param) + { + const auto atol = param.getDefault("atol", 1.0e-8); + const auto rtol = param.getDefault("rtol", 5.0e-12); + + return ErrorTolerance{ atol, rtol }; + } + + std::vector + loadReference(const ::Opm::parameter::ParameterGroup& param) + { + namespace fs = boost::filesystem; + + const auto fname = + fs::path(param.get("ref-dir")) / "trans.txt"; + + fs::ifstream input(fname); + + if (! input) { + std::ostringstream os; + + os << "Unable to Open Reference Trans-Data File " + << fname.filename(); + + throw std::domain_error(os.str()); + } + + return { + std::istream_iterator(input), + std::istream_iterator() + }; + } + + bool transfieldAcceptable(const ::Opm::parameter::ParameterGroup& param, + const std::vector& trans) + { + const auto Tref = loadReference(param); + + if (Tref.size() != trans.size()) { + return false; + } + + const auto diff = VectorDifference{ trans, Tref }; + + using Vector1 = std::decay::type; + using Vector2 = std::decay::type; + using Ratio = VectorRatio; + + const auto rat = Ratio(diff, Tref); + const auto tol = testTolerances(param); + + return ! ((pointMetric(diff) > tol.absolute) || + (pointMetric(rat) > tol.relative)); + } +} // namespace Anonymous + +int main(int argc, char* argv[]) +try { + const auto prm = example::initParam(argc, argv); + const auto pth = example::FilePaths(prm); + const auto G = example::initGraph(pth); + const auto T = G.transmissibility(); + const auto ok = transfieldAcceptable(prm, T); + + std::cout << (ok ? "OK" : "FAIL") << '\n'; + + if (! ok) { + return EXIT_FAILURE; + } +} +catch (const std::exception& e) { + std::cerr << "Caught Exception: " << e.what() << '\n'; + + return EXIT_FAILURE; +} diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclunithandling.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclunithandling.cpp new file mode 100644 index 0000000000..94bd19594d --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclunithandling.cpp @@ -0,0 +1,235 @@ +/* + Copyright 2017 SINTEF ICT, Applied Mathematics. + Copyright 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 // HAVE_CONFIG_H + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE + +#define BOOST_TEST_MODULE TEST_ASSEMBLED_CONNECTIONS + +#include +#include +#include + +#include + +#include +#include + +BOOST_AUTO_TEST_SUITE (Basic_Conversion) + +BOOST_AUTO_TEST_CASE (Constructor) +{ + auto M = ::Opm::ECLUnits::createUnitSystem(1); // METRIC + auto F = ::Opm::ECLUnits::createUnitSystem(2); // FIELD + auto L = ::Opm::ECLUnits::createUnitSystem(3); // LAB + auto P = ::Opm::ECLUnits::createUnitSystem(4); // PVT-M + + BOOST_CHECK_THROW(::Opm::ECLUnits::createUnitSystem( 5), std::runtime_error); + BOOST_CHECK_THROW(::Opm::ECLUnits::createUnitSystem(-1), std::runtime_error); +} + +BOOST_AUTO_TEST_CASE (Metric) +{ + auto M = ::Opm::ECLUnits::createUnitSystem(1); + + // Pressure (bars) + { + const auto scale = M->pressure(); + const auto expect = 100.0e3; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Reservoir Volume Rate (rm^3 / day) + { + const auto scale = M->reservoirRate(); + const auto expect = 1.157407407407407e-05; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Reservoir Volume (rm^3) + { + const auto scale = M->reservoirVolume(); + const auto expect = 1.0; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Time (day) + { + const auto scale = M->time(); + const auto expect = 86.400e+03; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Transmissibility ((cP * m^3) / (day * barsa)) + { + const auto scale = M->transmissibility(); + const auto expect = 1.157407407407407e-13; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE (Field) +{ + auto F = ::Opm::ECLUnits::createUnitSystem(2); + + // Pressure (psi) + { + const auto scale = F->pressure(); + const auto expect = 6.894757293168360e+03; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Reservoir Volume Rate (rb / day) + { + const auto scale = F->reservoirRate(); + const auto expect = 1.840130728333334e-06; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Reservoir Volume (rb) + { + const auto scale = F->reservoirVolume(); + const auto expect = 1.589872949280001e-01; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Time (day) + { + const auto scale = F->time(); + const auto expect = 86.400e+03; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Transmissibility ((cP * rb) / (day * psia)) + { + const auto scale = F->transmissibility(); + const auto expect = 2.668883979653090e-13; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE (Lab) +{ + auto L = ::Opm::ECLUnits::createUnitSystem(3); + + // Pressure (atm) + { + const auto scale = L->pressure(); + const auto expect = 101.325e+03; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Reservoir Volume Rate (r(cm)^3 / h) + { + const auto scale = L->reservoirRate(); + const auto expect = 2.777777777777778e-10; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Reservoir Volume (r(cm)^3) + { + const auto scale = L->reservoirVolume(); + const auto expect = 1.0e-06; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Time (hour) + { + const auto scale = L->time(); + const auto expect = 3600.0; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Transmissibility ((cP * (cm)^3) / (h * atm)) + { + const auto scale = L->transmissibility(); + const auto expect = 2.741453518655592e-18; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE (PVT_M) +{ + auto P = ::Opm::ECLUnits::createUnitSystem(4); + + // Pressure (atm) + { + const auto scale = P->pressure(); + const auto expect = 101.325e+03; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Reservoir Volume Rate (rm^3 / day) + { + const auto scale = P->reservoirRate(); + const auto expect = 1.157407407407407e-05; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Reservoir Volume (rm^3) + { + const auto scale = P->reservoirVolume(); + const auto expect = 1.0; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Time (day) + { + const auto scale = P->time(); + const auto expect = 86.400e+03; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } + + // Transmissibility ((cP * rm^3 / (day * atm)) + { + const auto scale = P->transmissibility(); + const auto expect = 1.142272299439830e-13; + + BOOST_CHECK_CLOSE(scale, expect, 1.0e-10); + } +} + +BOOST_AUTO_TEST_SUITE_END () diff --git a/ThirdParty/custom-opm-flowdiag-app/opmCore/opm/parser/eclipse/Units/Units.hpp b/ThirdParty/custom-opm-flowdiag-app/opmCore/opm/parser/eclipse/Units/Units.hpp index 33de0ed443..cc1dde7e66 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opmCore/opm/parser/eclipse/Units/Units.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opmCore/opm/parser/eclipse/Units/Units.hpp @@ -13,47 +13,60 @@ //=========================================================================== /* - Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010, 2011, 2012 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 . +Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. +Copyright 2009, 2010, 2011, 2012 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 . */ #ifndef OPM_UNITS_HEADER #define OPM_UNITS_HEADER - - /** - * \file - * Constants and routines to assist in handling units of measurement. These are - * geared towards handling common units in reservoir descriptions. - */ +The unit sets emplyed in ECLIPSE, in particular the FIELD units, +are quite inconsistent. Ideally one should choose units for a set +of base quantities like Mass,Time and Length and then derive the +units for e.g. pressure and flowrate in a consisten manner. However +that is not the case; for instance in the metric system we have: -namespace Opm -{ - namespace prefix +[Length] = meters +[time] = days +[mass] = kg + +This should give: + +[Pressure] = [mass] / ([length] * [time]^2) = kg / (m * days * days) + +Instead pressure is given in Bars. When it comes to FIELD units the +number of such examples is long. +*/ +namespace Opm { +namespace prefix /// Conversion prefix for units. - { - const double micro = 1.0e-6; /**< Unit prefix [\f$\mu\f$] */ - const double milli = 1.0e-3; /**< Unit prefix [m] */ - const double centi = 1.0e-2; /**< Non-standard unit prefix [c] */ - const double deci = 1.0e-1; /**< Non-standard unit prefix [d] */ - const double kilo = 1.0e3; /**< Unit prefix [k] */ - const double mega = 1.0e6; /**< Unit prefix [M] */ - const double giga = 1.0e9; /**< Unit prefix [G] */ - } // namespace prefix +{ +constexpr const double micro = 1.0e-6; /**< Unit prefix [\f$\mu\f$] */ +constexpr const double milli = 1.0e-3; /**< Unit prefix [m] */ +constexpr const double centi = 1.0e-2; /**< Non-standard unit prefix [c] */ +constexpr const double deci = 1.0e-1; /**< Non-standard unit prefix [d] */ +constexpr const double kilo = 1.0e3; /**< Unit prefix [k] */ +constexpr const double mega = 1.0e6; /**< Unit prefix [M] */ +constexpr const double giga = 1.0e9; /**< Unit prefix [G] */ +} // namespace prefix - namespace unit +namespace unit /// Definition of various units. /// All the units are defined in terms of international standard /// units (SI). Example of use: We define a variable \c k which @@ -68,162 +81,246 @@ namespace Opm /// using namespace Opm::prefix /// double k = 1.0*milli*darcy; /// \endcode - { - ///\name Common powers - /// @{ - inline double square(double v) { return v * v; } - inline double cubic (double v) { return v * v * v; } - /// @} +{ +///\name Common powers +/// @{ +constexpr double square(double v) { return v * v; } +constexpr double cubic (double v) { return v * v * v; } +/// @} - // -------------------------------------------------------------- - // Basic (fundamental) units and conversions - // -------------------------------------------------------------- +// -------------------------------------------------------------- +// Basic (fundamental) units and conversions +// -------------------------------------------------------------- - /// \name Length - /// @{ - const double meter = 1; - const double inch = 2.54 * prefix::centi*meter; - const double feet = 12 * inch; - /// @} +/// \name Length +/// @{ +constexpr const double meter = 1; +constexpr const double inch = 2.54 * prefix::centi*meter; +constexpr const double feet = 12 * inch; +/// @} - /// \name Time - /// @{ - const double second = 1; - const double minute = 60 * second; - const double hour = 60 * minute; - const double day = 24 * hour; - const double year = 365 * day; - /// @} +/// \name Time +/// @{ +constexpr const double second = 1; +constexpr const double minute = 60 * second; +constexpr const double hour = 60 * minute; +constexpr const double day = 24 * hour; +constexpr const double year = 365 * day; +/// @} - /// \name Volume - /// @{ - const double gallon = 231 * cubic(inch); - const double stb = 42 * gallon; - const double liter = 1 * cubic(prefix::deci*meter); - /// @} +/// \name Volume +/// @{ +constexpr const double gallon = 231 * cubic(inch); +constexpr const double stb = 42 * gallon; +constexpr const double liter = 1 * cubic(prefix::deci*meter); +/// @} - /// \name Mass - /// @{ - const double kilogram = 1; - // http://en.wikipedia.org/wiki/Pound_(mass)#Avoirdupois_pound - const double pound = 0.45359237 * kilogram; - /// @} +/// \name Mass +/// @{ +constexpr const double kilogram = 1; +constexpr const double gram = 1.0e-3 * kilogram; +// http://en.wikipedia.org/wiki/Pound_(mass)#Avoirdupois_pound +constexpr const double pound = 0.45359237 * kilogram; +/// @} - // -------------------------------------------------------------- - // Standardised constants - // -------------------------------------------------------------- +// -------------------------------------------------------------- +// Standardised constants +// -------------------------------------------------------------- - /// \name Standardised constant - /// @{ - const double gravity = 9.80665 * meter/square(second); - /// @} +/// \name Standardised constant +/// @{ +constexpr const double gravity = 9.80665 * meter/square(second); +/// @} - // -------------------------------------------------------------- - // Derived units and conversions - // -------------------------------------------------------------- +// -------------------------------------------------------------- +// Derived units and conversions +// -------------------------------------------------------------- - /// \name Force - /// @{ - const double Newton = kilogram*meter / square(second); // == 1 - const double lbf = pound * gravity; // Pound-force - /// @} +/// \name Force +/// @{ +constexpr const double Newton = kilogram*meter / square(second); // == 1 +constexpr const double dyne = 1e-5*Newton; +constexpr const double lbf = pound * gravity; // Pound-force + /// @} - /// \name Pressure - /// @{ - const double Pascal = Newton / square(meter); // == 1 - const double barsa = 100000 * Pascal; - const double atm = 101325 * Pascal; - const double psia = lbf / square(inch); - /// @} + /// \name Pressure + /// @{ +constexpr const double Pascal = Newton / square(meter); // == 1 +constexpr const double barsa = 100000 * Pascal; +constexpr const double atm = 101325 * Pascal; +constexpr const double psia = lbf / square(inch); +/// @} - /// \name Viscosity - /// @{ - const double Pas = Pascal * second; // == 1 - const double Poise = prefix::deci*Pas; - /// @} +/// \name Temperature. This one is more complicated +/// because the unit systems used by Eclipse (i.e. degrees +/// Celsius and degrees Fahrenheit require to add or +/// subtract an offset for the conversion between from/to +/// Kelvin +/// @{ +constexpr const double degCelsius = 1.0; // scaling factor °C -> K +constexpr const double degCelsiusOffset = 273.15; // offset for the °C -> K conversion - namespace perm_details { - const double p_grad = atm / (prefix::centi*meter); - const double area = square(prefix::centi*meter); - const double flux = cubic (prefix::centi*meter) / second; - const double velocity = flux / area; - const double visc = prefix::centi*Poise; - const double darcy = (velocity * visc) / p_grad; - // == 1e-7 [m^2] / 101325 - // == 9.869232667160130e-13 [m^2] - } - /// \name Permeability - /// @{ - /// - /// A porous medium with a permeability of 1 darcy permits a flow (flux) - /// of \f$1\,\mathit{cm}^3/s\f$ of a fluid with viscosity - /// \f$1\,\mathit{cP}\f$ (\f$1\,mPa\cdot s\f$) under a pressure gradient - /// of \f$1\,\mathit{atm}/\mathit{cm}\f$ acting across an area of - /// \f$1\,\mathit{cm}^2\f$. - /// - const double darcy = perm_details::darcy; - /// @} +constexpr const double degFahrenheit = 5.0/9; // scaling factor °F -> K +constexpr const double degFahrenheitOffset = 255.37; // offset for the °C -> K conversion + /// @} - /** - * Unit conversion routines. - */ - namespace convert { - /** - * Convert from external units of measurements to equivalent - * internal units of measurements. Note: The internal units of - * measurements are *ALWAYS*, and exclusively, SI. - * - * Example: Convert a double @c kx, containing a permeability value - * in units of milli-darcy (mD) to the equivalent value in SI units - * (i.e., \f$m^2\f$). - * \code - * using namespace Opm::unit; - * using namespace Opm::prefix; - * convert::from(kx, milli*darcy); - * \endcode - * - * @param[in] q Physical quantity. - * @param[in] unit Physical unit of measurement. - * @return Value of @c q in equivalent SI units of measurements. - */ - inline double from(const double q, const double unit) - { - return q * unit; - } + /// \name Viscosity + /// @{ +constexpr const double Pas = Pascal * second; // == 1 +constexpr const double Poise = prefix::deci*Pas; +/// @} - /** - * Convert from internal units of measurements to equivalent - * external units of measurements. Note: The internal units of - * measurements are *ALWAYS*, and exclusively, SI. - * - * Example: Convert a std::vector p, containing - * pressure values in the SI unit Pascal (i.e., unit::Pascal) to the - * equivalent values in Psi (unit::psia). - * \code - * using namespace Opm::unit; - * std::transform(p.begin(), p.end(), p.begin(), - * boost::bind(convert::to, _1, psia)); - * \endcode - * - * @param[in] q Physical quantity, measured in SI units. - * @param[in] unit Physical unit of measurement. - * @return Value of @c q in unit unit. - */ - inline double to(const double q, const double unit) - { - return q / unit; - } - } // namespace convert +namespace perm_details { +constexpr const double p_grad = atm / (prefix::centi*meter); +constexpr const double area = square(prefix::centi*meter); +constexpr const double flux = cubic (prefix::centi*meter) / second; +constexpr const double velocity = flux / area; +constexpr const double visc = prefix::centi*Poise; +constexpr const double darcy = (velocity * visc) / p_grad; +// == 1e-7 [m^2] / 101325 +// == 9.869232667160130e-13 [m^2] +} +/// \name Permeability +/// @{ +/// +/// A porous medium with a permeability of 1 darcy permits a flow (flux) +/// of \f$1\,\mathit{cm}^3/s\f$ of a fluid with viscosity +/// \f$1\,\mathit{cP}\f$ (\f$1\,mPa\cdot s\f$) under a pressure gradient +/// of \f$1\,\mathit{atm}/\mathit{cm}\f$ acting across an area of +/// \f$1\,\mathit{cm}^2\f$. +/// +constexpr const double darcy = perm_details::darcy; +/// @} + +/** +* Unit conversion routines. +*/ +namespace convert { +/** +* Convert from external units of measurements to equivalent +* internal units of measurements. Note: The internal units of +* measurements are *ALWAYS*, and exclusively, SI. +* +* Example: Convert a double @c kx, containing a permeability value +* in units of milli-darcy (mD) to the equivalent value in SI units +* (i.e., \f$m^2\f$). +* \code +* using namespace Opm::unit; +* using namespace Opm::prefix; +* convert::from(kx, milli*darcy); +* \endcode +* +* @param[in] q Physical quantity. +* @param[in] unit Physical unit of measurement. +* @return Value of @c q in equivalent SI units of measurements. +*/ +constexpr double from(const double q, const double unit) +{ + return q * unit; +} + +/** +* Convert from internal units of measurements to equivalent +* external units of measurements. Note: The internal units of +* measurements are *ALWAYS*, and exclusively, SI. +* +* Example: Convert a std::vector p, containing +* pressure values in the SI unit Pascal (i.e., unit::Pascal) to the +* equivalent values in Psi (unit::psia). +* \code +* using namespace Opm::unit; +* std::transform(p.begin(), p.end(), p.begin(), +* boost::bind(convert::to, _1, psia)); +* \endcode +* +* @param[in] q Physical quantity, measured in SI units. +* @param[in] unit Physical unit of measurement. +* @return Value of @c q in unit unit. +*/ +constexpr double to(const double q, const double unit) +{ + return q / unit; +} +} // namespace convert +} + +namespace Metric { +using namespace prefix; +using namespace unit; +constexpr const double Pressure = barsa; +constexpr const double Temperature = degCelsius; +constexpr const double TemperatureOffset = degCelsiusOffset; +constexpr const double AbsoluteTemperature = degCelsius; // actually [K], but the these two are identical +constexpr const double Length = meter; +constexpr const double Time = day; +constexpr const double Mass = kilogram; +constexpr const double Permeability = milli*darcy; +constexpr const double Transmissibility = centi*Poise*cubic(meter)/(day*barsa); +constexpr const double LiquidSurfaceVolume = cubic(meter); +constexpr const double GasSurfaceVolume = cubic(meter); +constexpr const double ReservoirVolume = cubic(meter); +constexpr const double GasDissolutionFactor = GasSurfaceVolume/LiquidSurfaceVolume; +constexpr const double OilDissolutionFactor = LiquidSurfaceVolume/GasSurfaceVolume; +constexpr const double Density = kilogram/cubic(meter); +constexpr const double PolymerDensity = kilogram/cubic(meter); +constexpr const double Salinity = kilogram/cubic(meter); +constexpr const double Viscosity = centi*Poise; +constexpr const double Timestep = day; +constexpr const double SurfaceTension = dyne/(centi*meter); +} -#ifndef HAS_ATTRIBUTE_UNUSED - namespace detail { - // Some units are sometimes unused, and generate a (potentially) large number of warnings - // Adding them here silences these warnings, and should have no side-effects - //JJS double __attribute__((unused)) unused_units = stb + liter + barsa + psia + darcy; - } // namespace detail -#endif +namespace Field { +using namespace prefix; +using namespace unit; +constexpr const double Pressure = psia; +constexpr const double Temperature = degFahrenheit; +constexpr const double TemperatureOffset = degFahrenheitOffset; +constexpr const double AbsoluteTemperature = degFahrenheit; // actually [°R], but the these two are identical +constexpr const double Length = feet; +constexpr const double Time = day; +constexpr const double Mass = pound; +constexpr const double Permeability = milli*darcy; +constexpr const double Transmissibility = centi*Poise*stb/(day*psia); +constexpr const double LiquidSurfaceVolume = stb; +constexpr const double GasSurfaceVolume = 1000*cubic(feet); +constexpr const double ReservoirVolume = stb; +constexpr const double GasDissolutionFactor = GasSurfaceVolume/LiquidSurfaceVolume; +constexpr const double OilDissolutionFactor = LiquidSurfaceVolume/GasSurfaceVolume; +constexpr const double Density = pound/cubic(feet); +constexpr const double PolymerDensity = pound/stb; +constexpr const double Salinity = pound/stb; +constexpr const double Viscosity = centi*Poise; +constexpr const double Timestep = day; +constexpr const double SurfaceTension = dyne/(centi*meter); +} - } // namespace unit -} // namespace Opm -#endif // OPM_UNITS_HEADER \ No newline at end of file + +namespace Lab { +using namespace prefix; +using namespace unit; +constexpr const double Pressure = atm; +constexpr const double Temperature = degCelsius; +constexpr const double TemperatureOffset = degCelsiusOffset; +constexpr const double AbsoluteTemperature = degCelsius; // actually [K], but the these two are identical +constexpr const double Length = centi*meter; +constexpr const double Time = hour; +constexpr const double Mass = gram; +constexpr const double Permeability = milli*darcy; +constexpr const double Transmissibility = centi*Poise*cubic(centi*meter)/(hour*atm); +constexpr const double LiquidSurfaceVolume = cubic(centi*meter); +constexpr const double GasSurfaceVolume = cubic(centi*meter); +constexpr const double ReservoirVolume = cubic(centi*meter); +constexpr const double GasDissolutionFactor = GasSurfaceVolume/LiquidSurfaceVolume; +constexpr const double OilDissolutionFactor = LiquidSurfaceVolume/GasSurfaceVolume; +constexpr const double Density = gram/cubic(centi*meter); +constexpr const double PolymerDensity = gram/cubic(centi*meter); +constexpr const double Salinity = gram/cubic(centi*meter); +constexpr const double Viscosity = centi*Poise; +constexpr const double Timestep = hour; +constexpr const double SurfaceTension = dyne/(centi*meter); +} + +} + +#endif // OPM_UNITS_HEADER 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 e99dd7ca4b..694738a293 100644 --- a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/TracerTofSolver.cpp +++ b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/TracerTofSolver.cpp @@ -323,18 +323,14 @@ namespace FlowDiagnostics // Compute time-of-flight and tracer. tracer_[cell] = is_start_[cell] ? 1.0 : upwind_tracer_contrib / total_influx; - //tof_[cell] = (pv_[cell]*tracer_[cell] + upwind_tof_contrib) / (total_influx*tracer_[cell]); - if ( tracer_[cell] > 0.0 ) - { + if (tracer_[cell] > 0.0) { tof_[cell] = (pv_[cell]*tracer_[cell] + upwind_tof_contrib) - / (total_influx * tracer_[cell]); + / (total_influx * tracer_[cell]); } - else - { + else { tof_[cell] = max_tof_; } - }