diff --git a/ApplicationCode/CMakeLists.txt b/ApplicationCode/CMakeLists.txt index 1290b1b1cf..3a6e3b141f 100644 --- a/ApplicationCode/CMakeLists.txt +++ b/ApplicationCode/CMakeLists.txt @@ -33,6 +33,8 @@ include_directories( ${ResInsight_SOURCE_DIR}/ThirdParty/custom-opm-common/opm-common/ ${custom-opm-parser_SOURCE_DIR}/opm-parser/ ${custom-opm-flowdiagnostics_SOURCE_DIR}/opm-flowdiagnostics/ + ${custom-opm-flowdiagnostics-applications_SOURCE_DIR}/opm-flowdiagnostics-applications/ + ${custom-opm-flowdiagnostics-applications_SOURCE_DIR}/opmCore ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/Adm @@ -284,6 +286,7 @@ endif() set( LINK_LIBRARIES custom-opm-parser custom-opm-flowdiagnostics + custom-opm-flowdiagnostics-applications WellPathImportSsihub diff --git a/ApplicationCode/UnitTests/opm-flowdiagnostics-Test.cpp b/ApplicationCode/UnitTests/opm-flowdiagnostics-Test.cpp index f2dd723b8c..7da675a1d5 100644 --- a/ApplicationCode/UnitTests/opm-flowdiagnostics-Test.cpp +++ b/ApplicationCode/UnitTests/opm-flowdiagnostics-Test.cpp @@ -2,9 +2,17 @@ #include #include +#include TEST(opm_flowdiagnostics_test, basic_construction) { auto g = Opm::AssembledConnections{}; auto s = Opm::FlowDiagnostics::CellSet{}; + try { + auto eg = Opm::ECLGraph::load("hei", "hopp"); + } + catch(const std::exception& e) + { + std::cerr << "Caught exception: " << e.what() << '\n'; + } } diff --git a/CMakeLists.txt b/CMakeLists.txt index e79c04e8bb..a14fb334a7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -208,7 +208,7 @@ endif(RESINSIGHT_ERT_EXTERNAL_LIB_ROOT OR RESINSIGHT_ERT_EXTERNAL_INCLUDE_ROOT) add_subdirectory(ThirdParty/custom-opm-parser) add_subdirectory(ThirdParty/custom-opm-parser/custom-opm-parser-tests) add_subdirectory(ThirdParty/custom-opm-flowdiagnostics) - +add_subdirectory(ThirdParty/custom-opm-flowdiagnostics-applications) ################################################################################ # NRLib diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/CMakeLists.txt b/ThirdParty/custom-opm-flowdiagnostics-applications/CMakeLists.txt new file mode 100644 index 0000000000..894140eeb5 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/CMakeLists.txt @@ -0,0 +1,27 @@ +cmake_minimum_required (VERSION 2.8) + +project (custom-opm-flowdiagnostics-applications) + +include_directories( + ../custom-opm-flowdiagnostics/opm-flowdiagnostics + opm-flowdiagnostics-applications + opmCore + ${ERT_INCLUDE_DIRS} + ${Boost_INCLUDE_DIRS} +) + +include (opm-flowdiagnostics-applications/CMakeLists_files.cmake) + +set(project_source_files + ${MAIN_SOURCE_FILES} + ${PUBLIC_HEADER_FILES} +) + +foreach (file ${project_source_files} ) + list(APPEND project_source_files_complete_path1 "opm-flowdiagnostics-applications/${file}" ) +endforeach() + +add_library(custom-opm-flowdiagnostics-applications + ${project_source_files_complete_path1} +) + diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/.gitignore b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/.gitignore new file mode 100644 index 0000000000..032094fc9f --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/.gitignore @@ -0,0 +1,66 @@ +# Compiled Object files +*.slo +*.lo +*.o +*.obj + +# Precompiled Headers +*.gch +*.pch + +# Compiled Dynamic libraries +*.so +*.dylib +*.dll + +# Fortran module files +*.mod + +# Compiled Static libraries +*.lai +*.la +*.a +*.lib + +# Executables +*.exe +*.out +*.app + +# Editor backup files +*~ +.\#* +\#*\# +.\#*\# + +# Eclipse project settings +.cproject +.project +.settings/* + +# QtCreator project settings +CMakeLists.txt.user* + +# In-tree build with CMake +CMakeCache.txt +CMakeFiles/ +cmake_install.cmake +config.h +opm-flowdiagnostics-applications-config.cmake +opm-flowdiagnostics-applications-config-version.cmake +opm-flowdiagnostics-applications-install.cmake +Makefile +bin/ +lib/ +Doxyfile +Documentation/html +*.pc +install_manifest.txt + +# Testing framework +CTestTestfile.cmake +DartConfiguration.tcl +Testing/ +.DS_Store +.idea +build diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/CMakeLists.txt b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/CMakeLists.txt new file mode 100644 index 0000000000..d6ca46957f --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/CMakeLists.txt @@ -0,0 +1,93 @@ +# -*- mode: cmake; tab-width: 2; indent-tabs-mode: t; truncate-lines: t; compile-command: "cmake -Wdev" -*- +# vim: set filetype=cmake autoindent tabstop=2 shiftwidth=2 noexpandtab softtabstop=2 nowrap: + +########################################################################### +# # +# Note: The bulk of the build system is located in the cmake/ directory. # +# This file only contains the specializations for this particular # +# project. Most likely you are interested in editing one of these # +# files instead: # +# # +# dune.module Name and version number # +# CMakeLists_files.cmake Path of source files # +# cmake/Modules/${project}-prereqs.cmake Dependencies # +# # +########################################################################### + +cmake_minimum_required (VERSION 2.8) + +# additional search modules +set(OPM_COMMON_ROOT "" CACHE PATH "Root directory containing OPM related cmake modules") +option(SIBLING_SEARCH "Search for other modules in sibling directories?" ON) + +if (NOT OPM_COMMON_ROOT) + find_package(opm-common QUIET) +endif() + +if (opm-common_FOUND) + include(OpmInit) +else() + unset(opm-common_FOUND) + + if (NOT OPM_COMMON_ROOT AND SIBLING_SEARCH) + set(OPM_COMMON_ROOT ${PROJECT_SOURCE_DIR}/../opm-common) + endif() + + if (OPM_COMMON_ROOT) + list(APPEND CMAKE_MODULE_PATH "${OPM_COMMON_ROOT}/cmake/Modules") + include (OpmInit OPTIONAL RESULT_VARIABLE OPM_INIT) + set(OPM_MACROS_ROOT ${OPM_COMMON_ROOT}) + endif() + + if (NOT OPM_INIT) + message("" ) + message(" /---------------------------------------------------------------------------------\\") + message(" | Could not locate the opm build macros. The opm build macros |") + message(" | are in a separate repository - instructions to proceed: |") + message(" | |") + message(" | 1. Clone the repository: git clone git@github.com:OPM/opm-common.git |") + message(" | |") + message(" | 2. Run cmake in the current project with -DOPM_COMMON_ROOT=/opm-common |") + message(" | |") + message(" \\---------------------------------------------------------------------------------/") + message("" ) + message(FATAL_ERROR "Could not find OPM Macros") + endif() +endif() + +# not the same location as most of the other projects; this hook overrides +macro (dir_hook) +endmacro (dir_hook) + +# 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 +include (${project}-prereqs) + +# read the list of components from this file (in the project directory); +# it should set various lists with the names of the files to include +include (CMakeLists_files.cmake) + +macro (config_hook) +endmacro (config_hook) + +macro (prereqs_hook) +endmacro (prereqs_hook) + +macro (sources_hook) +endmacro (sources_hook) + +macro (fortran_hook) +endmacro (fortran_hook) + +macro (files_hook) +endmacro (files_hook) + +macro (tests_hook) +endmacro (tests_hook) + +macro (install_hook) +endmacro (install_hook) + +# all setup common to the OPM library modules is done here +include (OpmLibMain) diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/CMakeLists_files.cmake b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/CMakeLists_files.cmake new file mode 100644 index 0000000000..1b429ff30e --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/CMakeLists_files.cmake @@ -0,0 +1,38 @@ +# This file sets up five lists: +# MAIN_SOURCE_FILES List of compilation units which will be included in +# the library. If it isn't on this list, it won't be +# part of the library. Please try to keep it sorted to +# maintain sanity. +# +# TEST_SOURCE_FILES List of programs that will be run as unit tests. +# +# TEST_DATA_FILES Files from the source three that should be made +# available in the corresponding location in the build +# tree in order to run tests there. +# +# EXAMPLE_SOURCE_FILES Other programs that will be compiled as part of the +# build, but which is not part of the library nor is +# run as tests. +# +# PUBLIC_HEADER_FILES List of public header files that should be +# distributed together with the library. The source +# files can of course include other files than these; +# you should only add to this list if the *user* of +# the library needs it. + +list (APPEND MAIN_SOURCE_FILES + opm/utility/ECLGraph.cpp + opm/utility/ECLWellSolution.cpp + ) + +list (APPEND TEST_SOURCE_FILES + ) + +list (APPEND EXAMPLE_SOURCE_FILES + examples/computeToFandTracers.cpp + ) + +list (APPEND PUBLIC_HEADER_FILES + opm/utility/ECLGraph.hpp + opm/utility/ECLWellSolution.hpp + ) diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/LICENSE b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/LICENSE new file mode 100644 index 0000000000..9cecc1d466 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/LICENSE @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + {one line to give the program's name and a brief idea of what it does.} + Copyright (C) {year} {name of author} + + This program 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. + + This program 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 this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + {project} Copyright (C) {year} {fullname} + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/README.md b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/README.md new file mode 100644 index 0000000000..d8782caf03 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/README.md @@ -0,0 +1,2 @@ +# opm-flowdiagnostics-applications +Stand-Alone Utilities for Developing and Testing Flow Diagnostics Compuational Kernels diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/dune.module b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/dune.module new file mode 100644 index 0000000000..1285983fed --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/dune.module @@ -0,0 +1,13 @@ +#################################################################### +# Dune module information file: This file gets parsed by dunecontrol +# and by the CMake build scripts. +#################################################################### + +Module: opm-flowdiagnostics-applications +Description: flow diagnostics applications and examples +Version: 2016.10-pre +Label: 2016.10-pre +Maintainer: bard.skaflestad@sintef.no +MaintainerName: Bård Skaflestad +Url: http://opm-project.org +Depends: opm-common opm-flowdiagnostics opm-core diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/examples/computeToFandTracers.cpp b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/examples/computeToFandTracers.cpp new file mode 100644 index 0000000000..fd757d0f25 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/examples/computeToFandTracers.cpp @@ -0,0 +1,206 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace { + bool isFile(const boost::filesystem::path& p) + { + namespace fs = boost::filesystem; + + auto is_regular_file = [](const fs::path& pth) + { + return fs::exists(pth) && fs::is_regular_file(pth); + }; + + return is_regular_file(p) + || (fs::is_symlink(p) && + is_regular_file(fs::read_symlink(p))); + } + + boost::filesystem::path + deriveFileName(boost::filesystem::path file, + const std::vector& extensions) + { + for (const auto& ext : extensions) { + file.replace_extension(ext); + + if (isFile(file)) { + return file; + } + } + + const auto prefix = file.parent_path() / file.stem(); + + std::ostringstream os; + + os << "Unable to derive valid filename from model prefix " + << prefix.generic_string(); + + throw std::invalid_argument(os.str()); + } + + Opm::FlowDiagnostics::ConnectionValues + extractFluxField(const Opm::ECLGraph& G, const int step) + { + 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 phas = ConnVals::PhaseID{0}; + + for (const auto& p : { Opm::BlackoilPhases::Aqua , + Opm::BlackoilPhases::Liquid , + Opm::BlackoilPhases::Vapour }) + { + const auto pflux = G.flux(p, step); + + if (! pflux.empty()) { + assert (pflux.size() == nconn.total); + + auto conn = ConnVals::ConnID{0}; + + for (const auto& v : pflux) { + flux(conn, phas) = v; + + conn.id += 1; + } + } + + phas.id += 1; + } + + return flux; + } + + Opm::FlowDiagnostics::Toolbox + initialiseFlowDiagnostics(const Opm::ECLGraph& G, + const std::vector& well_fluxes, + const int step) + { + const auto connGraph = Opm::FlowDiagnostics:: + ConnectivityGraph{ static_cast(G.numCells()), + G.neighbours() }; + + using FDT = Opm::FlowDiagnostics::Toolbox; + + auto fl = extractFluxField(G, step); + const size_t num_conn = fl.numConnections(); + const size_t num_phases = fl.numPhases(); + for (size_t conn = 0; conn < num_conn; ++conn) { + using Co = Opm::FlowDiagnostics::ConnectionValues::ConnID; + using Ph = Opm::FlowDiagnostics::ConnectionValues::PhaseID; + for (size_t phase = 0; phase < num_phases; ++phase) { + fl(Co{conn}, Ph{phase}) /= 86400; // HACK! converting to SI. + } + } + + Opm::FlowDiagnostics::CellSetValues inflow; + for (const auto& well : well_fluxes) { + for (const auto& completion : well.completions) { + const int grid_index = completion.grid_index; + const auto& ijk = completion.ijk; + const int cell_index = G.activeCell(ijk, grid_index); + inflow.addCellValue(cell_index, completion.reservoir_inflow_rate); + } + } + + // Create the Toolbox. + auto tool = FDT{ connGraph }; + tool.assignPoreVolume(G.poreVolume()); + tool.assignConnectionFlux(fl); + tool.assignInflowFlux(inflow); + + return tool; + } +} // Anonymous namespace + +// Syntax (typical): +// computeToFandTracers case= step= + +int main(int argc, char* argv[]) +try { + // 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); + + // 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); + Opm::ECLWellSolution wsol(restart); + auto well_fluxes = wsol.solution(step, graph.numGrids()); + auto fdTool = initialiseFlowDiagnostics(graph, well_fluxes, step); + + // Solve for time of flight. + using FDT = Opm::FlowDiagnostics::Toolbox; + std::vector start; + auto sol = fdTool.computeInjectionDiagnostics(start); + const auto& tof = sol.fd.timeOfFlight(); + + // Write it to standard out. + std::cout.precision(16); + for (double t : tof) { + std::cout << t << '\n'; + } +} +catch (const std::exception& e) { + std::cerr << "Caught exception: " << e.what() << '\n'; +} diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp new file mode 100644 index 0000000000..b863a15e2b --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp @@ -0,0 +1,1848 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#if HAVE_CONFIG_H +#include +#endif + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +/// \file +/// +/// Implementation of \c ECLGraph interface. + +namespace { + namespace ECL { + using GridPtr = ::ERT::ert_unique_ptr; + using FilePtr = ::ERT::ert_unique_ptr; + + /// Internalise on-disk representation of ECLIPSE grid. + /// + /// \param[in] grid Name or prefix of on-disk representation of + /// ECLIPSE grid. If using a prefix, the loader + /// will consider both .EGRID and .GRID versions of + /// the input. + /// + /// \return Internalised ERT Grid representation. + GridPtr loadCase(const boost::filesystem::path& grid); + + /// Internalise on-disk representation of ECL file. + /// + /// \param[in] file Name of ECLIPSE output file. + /// + /// \return Internalised ERT file contents. + FilePtr loadFile(const boost::filesystem::path& file); + + /// Retrieve total number of grids managed by model's main grid. + /// + /// \param[in] G Main grid obtained from loadCase(). + /// + /// \return Total number of grids in \p G. Equal to 1 + total + /// number of LGRs in model. + int numGrids(const ecl_grid_type* G); + + /// Access individual grid by numeric index. + /// + /// \param[in] G Main grid obtained from loadCase(). + /// + /// \param[in] gridID Numeric index of requested grid. Zero for the + /// main grid (i.e., \p G itself) or positive for one of the + /// LGRs. Must be strictly less than \code numGrids(G) \endcode. + /// + /// \return Pointer to ECL grid corresponding to numeric ID. + const ecl_grid_type* + getGrid(const ecl_grid_type* G, const int gridID); + + /// Extract Cartesian dimensions of an ECL grid. + /// + /// \param[in] G ERT grid instance corresponding to the model's main + /// grid or one of its LGRs. Typically obtained from function + /// getGrid(). + /// + /// \return Cartesian dimensions of \p G. Corresponds to number of + /// cells in each cardinal direction in 3D depositional space. + std::array + cartesianDimensions(const ecl_grid_type* G); + + /// Retrieve global pore-volume vector from INIT source. + /// + /// Specialised tool needed to determine the active cells. + /// + /// \param[in] G ERT Grid representation. + /// + /// \param[in] init ERT representation of INIT source. + /// + /// \return Vector of pore-volumes for all global cells of \p G. + std::vector + getPVolVector(const ecl_grid_type* G, + const ecl_file_type* init, + const int grid_ID = 0); + + /// Extract non-neighbouring connections from ECLIPSE model + /// + /// \param[in] G ERT Grid representation corresponding to model's + /// main grid obtained directly from loadCase(). + /// + /// \param[in] init ERT representation of INIT source. + /// + /// \return Model's non-neighbouring connections, including those + /// between main and local grids. + std::vector + loadNNC(const ecl_grid_type* G, + const ecl_file_type* init); + + class CartesianGridData + { + public: + /// Constructor. + /// + /// \param[in] G ERT grid structure corresponding either to the + /// model's main grid or, if applicable, one of its LGRs. + /// + /// \param[in] init Internalised ERT representation of result + /// set's INIT file. + /// + /// \param[in] gridID Numeric identifier of this grid. Zero for + /// main grid, positive for LGRs. + CartesianGridData(const ecl_grid_type* G, + const ecl_file_type* init, + const int gridID); + + /// Retrieve number of active cells in graph. + std::size_t numCells() const; + + /// Retrive number of connections in graph. + std::size_t numConnections() const; + + /// Retrive 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. + const std::vector& neighbours() const; + + /// Retrive 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. + const std::vector& activePoreVolume() const; + + /// Retrieve ID of active cell from global ID. + int activeCell(const std::size_t globalCell) const; + + /// Retrieve ID of active cell from (I,J,K) index tuple. + int activeCell(const int i, const int j, const int k) const; + + /// Predicate for whether or not a particular active cell is + /// further subdivided by an LGR. + /// + /// \param[in] cellID Index of particular active cell in this + /// grid. + /// + /// \return Whether or not cell identified by grid-local active + /// index \p cellID is further subdivided by an LGR. + bool isSubdivided(const int cellID) const; + + /// Retrieve values of result set vector for all global cells in + /// grid. + /// + /// Mostly for implementing connectionData(). + /// + /// \param[in] src ECLIPSE result set. + /// + /// \param[in] vector Name of result set vector. + /// + /// \return Numerical values of result set vector, relative to + /// global cell numbering of this grid. + std::vector + cellData(const ecl_file_type* src, + const std::string& vector) const; + + /// Retrieve values of result set vector for all Cartesian + /// connections in grid. + /// + /// \param[in] src ECLIPSE result set. + /// + /// \param[in] vector Name prefix of result set vector (e.g., + /// "FLROIL" for oil flux (flow-rate of oil)). + /// + /// \return Numerical values of result set vector attributed to + /// all of the grid's Cartesian connections. + std::vector + connectionData(const ecl_file_type* src, + const std::string& vector) const; + + private: + /// Facility for deriving Cartesian neighbourship in a grid + /// (main or LGR) and for mapping result set vectors to grid's + /// canonical (global) cells. + class CartesianCells + { + public: + /// Canonical directions of Cartesian neighbours. + enum class Direction { I, J, K }; + + /// Constructor + /// + /// \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(). + CartesianCells(const ecl_grid_type* G, + const std::vector& pvol); + + /// Retrive global cell indices of all active cells in grid. + std::vector activeGlobal() const; + + const std::vector& activePoreVolume() const; + + /// Map input vector to all global cells. + /// + /// \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 mapped to global cells or unchanged + /// if input is defined on some other subset. + template + std::vector + scatterToGlobal(const std::vector& x) const; + + /// Retrieve total number of cells in grid, including + /// inactive ones. + /// + /// Needed to allocate result vectors on global cells. + std::size_t numGlobalCells() const; + + /// Retrieve active cell ID of particular global cell. + /// + /// \param[in] globalCell Index of particular global cell. + /// + /// \return Active cell ID of \p globalCell. Returns + /// negative one (\code -1 \endcode) if \code globalCell >= + /// numGlobalCells \endcode or if the global cell is + /// inactive. + int getActiveCell(const std::size_t globalCell) const; + + /// Retrieve global cell ID of from (I,J,K) index tuple. + std::size_t + getGlobalCell(const int i, const int j, const int k) const; + + /// Retrieve active cell ID of particular global cell's + /// neighbour in given Cartesian direction. + /// + /// \param[in] globalCell Index of particular global cell. + /// + /// \param[in] d Cartesian direction in which to look for a + /// neighbouring cell. + /// + /// \return Active cell ID of \p globalCell's neighbour in + /// direction \d. Returns negative one (\code -1 \endcode) + /// if \code globalCell >= numGlobalCells \endcode or if the + /// global cell is inactive, or if there is no neighbour in + /// direction \p d (e.g., if purported neighbour would be + /// outside model). + int getNeighbour(const std::size_t globalCell, + const Direction d) const; + + /// Predicate for whether or not a particular active cell is + /// further subdivided by an LGR. + bool isSubdivided(const int cellID) const; + + private: + struct ResultSetMapping { + /// Explicit mapping between ACTNUM!=0 cells and global + /// cells. + struct ID { + std::size_t act; + std::size_t glob; + }; + + /// Number of explicitly active cells (SUM(ACTNUM != 0)). + std::size_t num_active; + + /// Active subset of global cells. + std::vector subset; + }; + + using IndexTuple = std::array; + + /// Size of grid's bounding box (i.e., the number of cells + /// in each cardinal direction in 3D depositional space). + const IndexTuple cartesianSize_; + + /// Map cell-based data vectors to grid's global cells. + ResultSetMapping rsMap_; + + /// Static pore-volumes of all active cells. + std::vector activePVol_; + + /// Active index of model's global cells. + std::vector active_ID_; + + /// 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); + + /// Compute linear index of global cell from explicit + /// (I,J,K) tuple. + /// + /// \param[in] ijk Explicit (I,J,K) tuple of global cell. + /// + /// \return Linear index (natural ordering) of global cell + /// (I,J,K). + std::size_t globIdx(const IndexTuple& ijk) const; + + /// Decompose global (linear) cell index into its (I,J,K) + /// index tuple. + /// + /// \param[in] globalCell Index of particular global cell. + /// Must be in the range \code [0 .. numGlobalCells()) + /// \endcode. + /// + /// \return Index triplet of \p globalCell's location within + /// model. + IndexTuple ind2sub(const std::size_t globalCell) const; + }; + + /// Collection of global (cell) IDs. + using GlobalIDColl = std::vector; + + /// Collection of (global) cell IDs corresponding to the flow + /// source of each connection. + using OutCell = + std::map; + + /// Collection of direction strings to simplify vector name + /// derivation (replaces chains of if-else) + using DirectionSuffix = + std::map; + + /// Numeric identity of this grid. Zero for main grid, greater + /// than zero for LGRs. + const int gridID_; + + /// Map results from active to global cells. + CartesianCells cells_; + + /// Known directional suffixes. + DirectionSuffix suffix_; + + /// Flattened neighbourship relation (array of size \code + /// 2*numConnections() \endcode). + std::vector neigh_; + + /// Source cells for each Cartesian connection. + OutCell outCell_; + + /// Predicate for whether or not a particular result vector is + /// defined on the grid's cells. + /// + /// \param[in] src Result set. + /// + /// \param[in] vector Name of result vector. + /// + /// \return Whether or not \p vector is defined on model's + /// cells and part of the result set \p src. + bool haveCellData(const ecl_file_type* src, + const std::string& vector) const; + + /// Predicate for whether or not a particular result vector is + /// defined on the grid's Cartesian connections. + /// + /// \param[in] src Result set. + /// + /// \param[in] vector Prefix of result vector name. + /// + /// \return Whether or not all vectors formed by \p vector plus + /// known directional suffixes are defined on model's cells and + /// part of the result set \p src. + bool haveConnData(const ecl_file_type* src, + const std::string& vector) const; + + /// Append directional cell data to global collection of + /// connection data identified by vector name prefix. + /// + /// \param[in] src Result set. + /// + /// \param[in] d Cartesian direction. + /// + /// \param[in] vector Prefix of result vector name. + /// + /// \param[in,out] x Global collection of connection data. On + /// input, collection of values corresponding to any previous + /// directions (preserved), and on output additionally contains + /// the data corresponding to the Cartesian direction \p d. + void connectionData(const ecl_file_type* src, + const CartesianCells::Direction d, + const std::string& vector, + std::vector& x) const; + + /// Form complete name of directional result set vector from + /// prefix and identified direction. + /// + /// \param[in] vector Prefix of result vector name. + /// + /// \param[in] d Cartesian direction. + /// + /// \return \code vector + suffix_[d] \endcode. + std::string + vectorName(const std::string& vector, + const CartesianCells::Direction d) const; + + /// Derive neighbourship relations on active cells in particular + /// Cartesian directions. + /// + /// Writes to \c neigh_ and \c outCell_. + /// + /// \param[in] gcells Collection of global (relative to \c + /// gridID_) cells that should be considered active (strictly + /// positive pore-volume and not deactivated through + /// ACTNUM=0). + /// + /// \param[in] init Internalised + void deriveNeighbours(const std::vector& gcells, + const ecl_file_type* init, + const CartesianCells::Direction d); + }; + } // namespace ECL +} // Anonymous namespace + +// ====================================================================== + +int ECL::numGrids(const ecl_grid_type* G) +{ + return 1 + ecl_grid_get_num_lgr(G); // Main + #LGRs. +} + +const ecl_grid_type* +ECL::getGrid(const ecl_grid_type* G, const int gridID) +{ + assert ((gridID >= 0) && "Grid ID must be non-negative"); + + if (gridID == ECL_GRID_MAINGRID_LGR_NR) { + return G; + } + + return ecl_grid_iget_lgr(G, gridID - 1); +} + +std::vector +ECL::getPVolVector(const ecl_grid_type* G, + const ecl_file_type* init, + const int gridID) +{ + auto make_szt = [](const int i) + { + return static_cast::size_type>(i); + }; + + const auto nglob = make_szt(ecl_grid_get_global_size(G)); + + auto pvol = std::vector(nglob, 1.0); + + if (ecl_file_has_kw(init, "PORV")) { + auto porv = + ecl_file_iget_named_kw(init, "PORV", gridID); + + assert ((make_szt(ecl_kw_get_size(porv)) == nglob) + && "Pore-volume must be provided for all global cells"); + + ecl_kw_get_data_as_double(porv, pvol.data()); + } + + return pvol; +} + +ECL::GridPtr +ECL::loadCase(const boost::filesystem::path& grid) +{ + auto G = GridPtr{ + ecl_grid_load_case(grid.generic_string().c_str()) + }; + + if (! G) { + std::ostringstream os; + + os << "Failed to load ECL Grid from '" + << grid.generic_string() << '\''; + + throw std::invalid_argument(os.str()); + } + + return G; +} + +ECL::FilePtr +ECL::loadFile(const boost::filesystem::path& file) +{ + // Read-only, keep open between requests + const auto open_flags = 0; + + auto F = FilePtr{ + ecl_file_open(file.generic_string().c_str(), open_flags) + }; + + if (! F) { + std::ostringstream os; + + os << "Failed to load ECL File object from '" + << file.generic_string() << '\''; + + throw std::invalid_argument(os.str()); + } + + return F; +} + +std::array +ECL::cartesianDimensions(const ecl_grid_type* G) +{ + auto make_szt = [](const int i) + { + return static_cast(i); + }; + + return { { make_szt(ecl_grid_get_nx(G)) , + make_szt(ecl_grid_get_ny(G)) , + make_szt(ecl_grid_get_nz(G)) } }; +} + +std::vector +ECL::loadNNC(const ecl_grid_type* G, + const ecl_file_type* init) +{ + auto make_szt = [](const int n) + { + return static_cast::size_type>(n); + }; + + auto nncData = std::vector{}; + + const auto numNNC = make_szt(ecl_nnc_export_get_size(G)); + + if (numNNC > 0) { + nncData.resize(numNNC); + + ecl_nnc_export(G, init, nncData.data()); + } + + return nncData; +} + +// ====================================================================== + +ECL::CartesianGridData:: +CartesianCells::CartesianCells(const ecl_grid_type* G, + const std::vector& pvol) + : cartesianSize_(::ECL::cartesianDimensions(G)) +{ + if (pvol.size() != static_cast + (this->cartesianSize_[0] * + this->cartesianSize_[1] * + this->cartesianSize_[2])) + { + throw std::invalid_argument("Grid must have PORV for all cells"); + } + + auto make_szt = [](const int i) + { + return static_cast(i); + }; + + using ID = ResultSetMapping::ID; + + this->rsMap_.num_active = make_szt(ecl_grid_get_nactive(G)); + + this->rsMap_.subset.clear(); + this->rsMap_.subset.reserve(this->rsMap_.num_active); + + for (decltype(ecl_grid_get_nactive(G)) + act = 0, nact = ecl_grid_get_nactive(G); + act < nact; ++act) + { + const auto glob = + make_szt(ecl_grid_get_global_index1A(G, act)); + + if (pvol[glob] > 0.0) { + this->rsMap_.subset.push_back(ID{ make_szt(act), glob }); + } + } + + { + std::vector(pvol.size(), -1).swap(this->active_ID_); + + this->activePVol_.clear(); + this->activePVol_.reserve(this->rsMap_.subset.size()); + + this->is_divided_.clear(); + this->is_divided_.reserve(this->rsMap_.subset.size()); + + auto active = 0; + + for (const auto& cell : this->rsMap_.subset) { + this->active_ID_[cell.glob] = active++; + this->activePVol_.push_back(pvol[cell.glob]); + + const auto ert_active = static_cast(cell.act); + const auto is_divided = + nullptr != ecl_grid_get_cell_lgr1A(G, ert_active); + + this->is_divided_.push_back(is_divided); + } + } +} + +std::vector +ECL::CartesianGridData::CartesianCells::activeGlobal() const +{ + auto active = std::vector{}; + active.reserve(this->rsMap_.subset.size()); + + for (const auto& id : this->rsMap_.subset) { + active.push_back(id.glob); + } + + return active; +} + +const std::vector& +ECL::CartesianGridData::CartesianCells::activePoreVolume() const +{ + return this->activePVol_; +} + +template +std::vector +ECL::CartesianGridData:: +CartesianCells::scatterToGlobal(const std::vector& x) const +{ + // Assume that input vector 'x' is either defined on explicit notion of + // active cells (ACTNUM != 0) or on all global cells or some other + // contiguous index set (e.g., the NNCs). + + const auto num_explicit_active = + static_cast(this->rsMap_.num_active); + + if (x.size() != num_explicit_active) { + // Input not defined on explictly active cells. Let caller deal + // with it. This typically corresponds to the set of global cells + // or the list of NNCs. + return x; + } + + auto y = std::vector(this->numGlobalCells()); + + for (const auto& i : this->rsMap_.subset) { + y[i.glob] = x[i.act]; + } + + return y; +} + +std::size_t +ECL::CartesianGridData::CartesianCells::numGlobalCells() const +{ + return this->active_ID_.size(); +} + +int +ECL::CartesianGridData:: +CartesianCells::getActiveCell(const std::size_t globalCell) const +{ + if (globalCell >= numGlobalCells()) { return -1; } + + return this->active_ID_[globalCell]; +} + +std::size_t +ECL::CartesianGridData:: +CartesianCells::getGlobalCell(const int i, const int j, const int k) const +{ + const auto ijk = IndexTuple { + static_cast(i), + static_cast(j), + static_cast(k), + }; + + return this->globIdx(ijk); +} + +int +ECL::CartesianGridData:: +CartesianCells::getNeighbour(const std::size_t globalCell, + const Direction d) const +{ + if (globalCell >= numGlobalCells()) { return -1; } + + auto ijk = ind2sub(globalCell); + + if (d == Direction::I) { ijk[0] += 1; } + else if (d == Direction::J) { ijk[1] += 1; } + else if (d == Direction::K) { ijk[2] += 1; } + else { + return -1; + } + + const auto globNeigh = globIdx(ijk); + + if (globNeigh >= numGlobalCells()) { return -1; } + + return this->active_ID_[globNeigh]; +} + +bool +ECL::CartesianGridData::CartesianCells::isSubdivided(const int cellID) const +{ + const auto ix = + static_castis_divided_.size())>(cellID); + + assert ((cellID >= 0) && (ix < this->is_divided_.size())); + + return this->is_divided_[ix]; +} + +std::size_t +ECL::CartesianGridData:: +CartesianCells::globIdx(const IndexTuple& ijk) const +{ + const auto& dim = this->cartesianSize_; + + for (auto d = 0*dim.size(), nd = dim.size(); d < nd; ++d) { + if (ijk[d] >= dim[d]) { return -1; } + } + + return ijk[0] + dim[0]*(ijk[1] + dim[1]*ijk[2]); +} + +ECL::CartesianGridData::CartesianCells::IndexTuple +ECL::CartesianGridData:: +CartesianCells::ind2sub(const std::size_t globalCell) const +{ + assert (globalCell < numGlobalCells()); + + auto ijk = IndexTuple{}; + auto g = globalCell; + + const auto& dim = this->cartesianSize_; + + ijk[0] = g % dim[0]; g /= dim[0]; + ijk[1] = g % dim[1]; + ijk[2] = g / dim[1]; assert (ijk[2] < dim[2]); + + assert (globIdx(ijk) == globalCell); + + return ijk; +} + +// ====================================================================== + +ECL::CartesianGridData::CartesianGridData(const ecl_grid_type* G, + const ecl_file_type* init, + const int gridID) + : gridID_(gridID) + , cells_ (G, ::ECL::getPVolVector(G, init, gridID_)) +{ + { + using VT = DirectionSuffix::value_type; + + suffix_.insert(VT(CartesianCells::Direction::I, "I+")); + suffix_.insert(VT(CartesianCells::Direction::J, "J+")); + suffix_.insert(VT(CartesianCells::Direction::K, "K+")); + } + + const auto gcells = this->cells_.activeGlobal(); + + // Too large, but this is a quick estimate. + this->neigh_.reserve(3 * (2 * this->numCells())); + + for (const auto d : { CartesianCells::Direction::I , + CartesianCells::Direction::J , + CartesianCells::Direction::K }) + { + this->deriveNeighbours(gcells, init, d); + } +} + +std::size_t +ECL::CartesianGridData::numCells() const +{ + return this->activePoreVolume().size(); +} + +std::size_t +ECL::CartesianGridData::numConnections() const +{ + return this->neigh_.size() / 2; +} + +const std::vector& +ECL::CartesianGridData::neighbours() const +{ + return this->neigh_; +} + +const std::vector& +ECL::CartesianGridData::activePoreVolume() const +{ + return this->cells_.activePoreVolume(); +} + +int +ECL::CartesianGridData::activeCell(const std::size_t globalCell) const +{ + return this->cells_.getActiveCell(globalCell); +} + +int +ECL::CartesianGridData::activeCell(const int i, + const int j, + const int k) const +{ + return this->activeCell(this->cells_.getGlobalCell(i, j, k)); +} + +bool +ECL::CartesianGridData::isSubdivided(const int cellID) const +{ + return this->cells_.isSubdivided(cellID); +} + +std::vector +ECL::CartesianGridData::cellData(const ecl_file_type* src, + const std::string& vector) const +{ + if (! this->haveCellData(src, vector)) { + return {}; + } + + const auto v = + ecl_file_iget_named_kw(src, vector.c_str(), this->gridID_); + + auto x = std::vector(ecl_kw_get_size(v)); + + ecl_kw_get_data_as_double(v, x.data()); + + return this->cells_.scatterToGlobal(x); +} + +bool +ECL::CartesianGridData::haveCellData(const ecl_file_type* src, + const std::string& vector) const +{ + // Recall: get_num_named_kw() is block aware (uses src->active_map). + + return ecl_file_get_num_named_kw(src, vector.c_str()) > this->gridID_; +} + +bool +ECL::CartesianGridData::haveConnData(const ecl_file_type* src, + const std::string& vector) const +{ + auto have_data = true; + + for (const auto& d : { CartesianCells::Direction::I , + CartesianCells::Direction::J , + CartesianCells::Direction::K }) + { + const auto vname = this->vectorName(vector, d); + have_data = this->haveCellData(src, vname); + + if (! have_data) { break; } + } + + return have_data; +} + +std::vector +ECL::CartesianGridData::connectionData(const ecl_file_type* src, + const std::string& vector) const +{ + if (! this->haveConnData(src, vector)) { + return {}; + } + + auto x = std::vector{}; x.reserve(this->numConnections()); + + for (const auto& d : { CartesianCells::Direction::I , + CartesianCells::Direction::J , + CartesianCells::Direction::K }) + { + this->connectionData(src, d, vector, x); + } + + return x; +} + +void +ECL::CartesianGridData:: +connectionData(const ecl_file_type* src, + const CartesianCells::Direction d, + const std::string& vector, + std::vector& x) const +{ + const auto v = this->cellData(src, this->vectorName(vector, d)); + + const auto& cells = this->outCell_.find(d); + + assert ((cells != this->outCell_.end()) && + "Direction must be I, J, or K"); + + for (const auto& cell : cells->second) { + x.push_back(v[cell]); + } +} + +std::string +ECL::CartesianGridData:: +vectorName(const std::string& vector, + const CartesianCells::Direction d) const +{ + const auto i = this->suffix_.find(d); + + assert ((i != this->suffix_.end()) && + "Direction must be I, J, or K"); + + return vector + i->second; +} + +void +ECL::CartesianGridData:: +deriveNeighbours(const std::vector& gcells, + const ecl_file_type* init, + const CartesianCells::Direction d) +{ + auto tran = std::string{"TRAN"}; + + switch (d) { + case CartesianCells::Direction::I: + tran += 'X'; + break; + + case CartesianCells::Direction::J: + tran += 'Y'; + break; + + case CartesianCells::Direction::K: + tran += 'Z'; + break; + + default: + throw std::invalid_argument("Input direction must be (I,J,K)"); + } + + const auto& T = this->haveCellData(init, tran) + ? this->cellData(init, tran) + : std::vector(this->cells_.numGlobalCells(), 1.0); + + auto& ocell = this->outCell_[d]; + ocell.reserve(gcells.size()); + + for (const auto& globID : gcells) { + const auto c1 = this->cells_.getActiveCell(globID); + + assert ((c1 >= 0) && "Internal error in active cell derivation"); + + if (this->cells_.isSubdivided(c1)) { + // Don't form connections to subdivided cells. We care only + // about the final refinement level (i.e., the most nested LGR + // object) and the connections are handled by the NNC code. + continue; + } + + if (T[globID] > 0.0) { + const auto other = this->cells_.getNeighbour(globID, d); + + if ((other >= 0) && ! this->cells_.isSubdivided(other)) { + assert (c1 != other); + + this->neigh_.push_back(c1); + this->neigh_.push_back(other); + + ocell.push_back(globID); + } + } + } +} + +// ===================================================================== + +/// Implementation of ECLGraph interface. +class Opm::ECLGraph::Impl +{ +public: + /// Constructor + /// + /// \param[in] grid Name or prefix of ECL grid (i.e., .GRID or + /// .EGRID) file. + /// + /// \param[in] init Name of ECL INIT file corresponding to \p grid + /// input. Assumed to provide at least a complete set + /// of pore-volume values (i.e., for all global cells + /// defined in the \p grid). + /// + /// If available in the INIT file, the constructor will + /// also leverage the transmissibility data when + /// constructing the active cell neighbourship table. + Impl(const Path& grid, const Path& init); + + /// Assign source object for phase flux calculation. + /// + /// \param[in] src Name of ECL restart file, possibly unified, from + /// which next set of phase fluxes should be retrieved. + void assignDataSource(const Path& src); + + /// Retrieve number of grids. + /// + /// \return The number of LGR grids plus one (the main grid). + int numGrids() const; + + /// Retrieve active cell ID from (I,J,K) tuple in particular grid. + /// + /// \param[in] gridID Identity of specific grid to which to relate the + /// (I,J,K) tuple. Zero for main grid and positive indices for any + /// LGRs. The (I,J,K) indices must be within the ranges implied by + /// the specific grid. + /// + /// \param[in] ijk Cartesian index tuple of particular cell. + /// + /// \return Active ID (relative to linear, global numbering) of cell \p + /// ijk from specified grid. Negative one (-1) if (I,J,K) outside + /// valid range or if the specific cell identified by \p ijk and \p + /// gridID is not actually active. + int activeCell(const int gridID, + const std::array& ijk) const; + + /// Retrieve number of active cells in graph. + std::size_t numCells() const; + + /// Retrive number of connections in graph. + std::size_t numConnections() const; + + /// Retrive 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. + /// + /// 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; + + /// Retrive phase flux on all connections defined by \code neighbours() + /// \endcode. + /// + /// \param[in] phase Canonical phase for which to retrive flux. + /// + /// \param[in] rptstep Selected temporal vector. Report-step ID. + /// + /// \return Flux values corresponding to selected phase and report step. + /// Empty if unavailable in the result set (e.g., by querying the gas + /// flux in an oil/water system or if the specified \p occurrence is not + /// reported due to report frequencies or no flux values are output at + /// all). + std::vector + flux(const BlackoilPhases::PhaseIndex phase, + const int rptstep) 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). + class NonNeighKeywordIndexSet + { + public: + /// Establish mapping between particular non-Cartesian neighbourship + /// relation and particular entry within a grid's keyword data. + struct Map { + /// Non-Cartesian neighbourship relation. + std::size_t neighIdx; + + /// Index into grid's keyword data. + std::size_t kwIdx; + }; + + using MapCollection = std::vector; + + /// Record a mapping in a particular grid. + /// + /// \param[in] grid Particular model grid for which to record a + /// mapping. + /// + /// \param[in] entry Individual index map. + void add(const int grid, Map&& entry); + + /// Retrieve collection of index maps for particular grid. + /// + /// \param[in] grid Specific model grid ID. Must be non-negative + /// and strictly less than the total number of grids in the + /// model. + /// + /// \return Collection of index maps attributable to \p grid. Empty + /// if no such collection exists. + const MapCollection& getGridCollection(const int grid) const; + + private: + using KWEntries = std::map; + + /// Collection of all index maps attributable to all grids for this + /// set of ECL keywords. + KWEntries subset_; + + /// Return value for the case of no existing collection. + MapCollection empty_; + }; + + /// Collection of all non-Cartesian neighbourship relations classified + /// according to connection type. + class NNC + { + public: + /// Classification of non-Cartesian neighbourship relations. + enum class Category { + /// Traditional non-neighbouring connections entirely internal + /// to a grid. Typically due to faults or fully unstructured + /// grid descriptions. Keywords NNC{1,2}, TRANNNC, and FLR*N+. + /// Positive from NNC1 to NNC2. + Normal, + + /// Connections between main grid and LGRs. Keywords NNCG, + /// NNCL, TRANGL, and FLR*L+. Positive from NNCG to NNCL. + GlobalToLocal, + + /// Connections between LGRs. Either due to two LGRs being + /// neighbouring entities in physical space or one LGR being + /// nested within another. Keywords NNA{1,2}, TRANLL, and + /// FLR*A+. Positive from NNA1 to NNA2. + Amalgamated, + }; + + /// Map a collection of non-Cartesian neighbourship relations to a + /// specific flux vector identifier. + struct FluxRelation { + /// Flux vector identifier. Should be one of "N+" for Normal + /// connections, "L+" for GlobalToLocal connections, and "A+" + /// for Amalgamated connections. + std::string fluxID; + + /// Collection of non-Cartesian neighbourship relations. + NonNeighKeywordIndexSet indexSet; + }; + + /// Constructor. + NNC(); + + /// Potentially record a new non-Cartesian connection. + /// + /// Will classify the connection according to the grids involved and + /// actually record the connection if both cells are active and + /// neither are subdivided. + /// + /// \param[in] grids Collection of all active grids in model. + /// + /// \param[in] offset Start index into global linear number for all + /// active grids. + /// + /// \param[in] nnc Non-neighbouring connection from result set. + void add(const std::vector& grids, + const std::vector& offset, + const ecl_nnc_type& nnc); + + std::vector allCategories() const; + + /// Retrieve total number of active non-neighbouring connections. + std::size_t numConnections() const; + + /// Access all active non-neighbouring connections. + const std::vector& getNeighbours() const; + + /// Retrieve all non-neighbouring connections of a particular + /// category (i.e., pertaining to a particular set of keywords). + /// + /// \param[in] type Category of non-neighbouring connections. + /// + /// \return All non-neighbouring connections of category \p type. + const FluxRelation& getRelations(const Category& type) const; + + private: + using KeywordIndexMap = std::map; + + /// Active non-Cartesian (non-neighbouring) connections. Cell IDs + /// in linear numbering of all model's active cells. + std::vector neigh_; + + /// Collection of + KeywordIndexMap keywords_; + + /// Factory for FluxRelations. + /// + /// Simplifies implementation of ctor. + /// + /// \param[in] cat Class + FluxRelation makeRelation(const Category cat) const; + + /// Identify connection category from connection's grids. + /// + /// - Normal connection if both grid IDs equal + /// - GlobalToLocal if one grid is the main model grid and the + /// other is an LGR. + /// - Amalgamated if both grids are LGRs. + /// + /// \param[in] grid1 Numeric identity of connection's first grid. + /// Zero if main grid, positive if LGR. + /// + /// \param[in] grid2 Numeric identity of connection's second grid. + /// Zero if main grid, positive if LGR. + /// + /// \return Category of connection between \p grid1 and \p grid2. + Category classifyConnection(const int grid1, const int grid2) const; + + /// Check if cell is viable connection endpoint in grid. + /// + /// A cell is a viable connection endpoint if it is active within a + /// particular grid and not further subdivided by an LGR. + /// + /// \param[in] grids Collection of all active grids in model. + /// + /// \param[in] gridID Numeric identity of connection grid. + /// Zero if main grid, positive if LGR. + /// + /// \param[in] cellID Global ID (relative to \p grid) of candidate + /// connection endpoint. + /// + /// \return Whether or not \p cellID is a viable connection endpoint + /// within \p gridID. + bool isViable(const std::vector& grids, + const int gridID, + const std::size_t cellID) const; + + /// Check if connection is viable + /// + /// A candidate non-Cartesian connection is viable if both of its + /// endpoints satisfy the viability criterion. + /// + /// \param[in] nnc Candidate non-Cartesian connection. + /// + /// \return Whether or not both candidate endpoints satisfy the + /// viability criterion. + bool isViable(const std::vector& grids, + const ecl_nnc_type& nnc) const; + }; + + /// Collection of model's non-neighbouring connections--be they within a + /// grid or between grids. + NNC nnc_; + + /// Collection of model's grids (main + LGRs). + std::vector grid_; + + /// Map each grid's active cellIDs to global numbering (in the index + /// range \code [0 .. numCells()) \endcode). + std::vector activeOffset_; + + /// Current result set. + ECL::FilePtr src_; + + /// Extract explicit non-neighbouring connections from ECL output. + /// + /// Writes to \c neigh_ and \c nncID_. + /// + /// \param[in] G ERT Grid representation. + /// + /// \param[in] init ERT representation of INIT source. + /// + /// \param[in] coll Backing data for neighbourship extraction. + void defineNNCs(const ecl_grid_type* G, + const ecl_file_type* init); + + /// Compute ECL vector basename for particular phase flux. + /// + /// \param[in] phase Canonical phase for which to derive ECL vector + /// basename. + /// + /// \return Basename for ECl vector corresponding to particular phase + /// flux. + std::string + flowVector(const BlackoilPhases::PhaseIndex phase) const; + + /// Extract flux values corresponding to particular result set vector + /// for all identified non-neighbouring connections. + /// + /// \param[in] vector Result set vector prefix. Typically computed by + /// method flowVector(). + /// + /// \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). + void fluxNNC(const std::string& vector, + std::vector& flux) const; +}; + +// ====================================================================== + +void +Opm::ECLGraph::Impl::NonNeighKeywordIndexSet:: +add(const int grid, Map&& entry) +{ + this->subset_[grid].push_back(std::move(entry)); +} + +const +Opm::ECLGraph::Impl::NonNeighKeywordIndexSet::MapCollection& +Opm::ECLGraph::Impl::NonNeighKeywordIndexSet:: +getGridCollection(const int grid) const +{ + auto coll = this->subset_.find(grid); + + if (coll == this->subset_.end()) { + // No NNCs of this category for this grid. Return empty. + return this->empty_; + } + + return coll->second; +} + +// ====================================================================== + +Opm::ECLGraph::Impl::NNC::NNC() +{ + using VT = KeywordIndexMap::value_type; + + for (const auto& cat : this->allCategories()) { + this->keywords_.insert(VT(cat, this->makeRelation(cat))); + } +} + +std::vector +Opm::ECLGraph::Impl::NNC::allCategories() const +{ + return { Category::Normal , + Category::GlobalToLocal , + Category::Amalgamated }; +} + +void +Opm::ECLGraph::Impl:: +NNC::add(const std::vector& grid, + const std::vector& offset, + const ecl_nnc_type& nnc) +{ + if (! this->isViable(grid, nnc)) { + // Zero transmissibility or at least one endpoint unviable. Don't + // record connection. + return; + } + + const auto neighIdx = this->numConnections(); + + { + const auto c = grid[nnc.grid_nr1].activeCell(nnc.global_index1); + const auto o = static_cast(offset[nnc.grid_nr1]); + + this->neigh_.push_back(o + c); + } + + { + const auto c = grid[nnc.grid_nr2].activeCell(nnc.global_index2); + const auto o = static_cast(offset[nnc.grid_nr2]); + + this->neigh_.push_back(o + c); + } + + const auto cat = this->classifyConnection(nnc.grid_nr1, nnc.grid_nr2); + + auto entry = NonNeighKeywordIndexSet::Map { + neighIdx, + static_cast(nnc.input_index) + }; + + this->keywords_[cat].indexSet.add(nnc.grid_nr2, std::move(entry)); +} + +std::size_t +Opm::ECLGraph::Impl::NNC::numConnections() const +{ + assert (this->neigh_.size() % 2 == 0); + + return this->neigh_.size() / 2; +} + +const std::vector& +Opm::ECLGraph::Impl::NNC::getNeighbours() const +{ + return this->neigh_; +} + +const Opm::ECLGraph::Impl::NNC::FluxRelation& +Opm::ECLGraph::Impl::NNC::getRelations(const Category& cat) const +{ + auto r = this->keywords_.find(cat); + + assert ((r != this->keywords_.end()) && + "Input category must be Normal, " + "GlobalToLocal or Amalgamated"); + + return r->second; +} + +Opm::ECLGraph::Impl::NNC::FluxRelation +Opm::ECLGraph::Impl::NNC::makeRelation(const Category cat) const +{ + switch (cat) { + case Category::Normal: + return { "N+", {} }; + + case Category::GlobalToLocal: + return { "L+", {} }; + + case Category::Amalgamated: + return { "A+", {} }; + } + + throw std::invalid_argument("Category must be Normal, " + "GlobalToLocal, or Amalgamated"); +} + +Opm::ECLGraph::Impl::NNC::Category +Opm::ECLGraph::Impl::NNC:: +classifyConnection(const int grid1, const int grid2) const +{ + if (grid1 == grid2) { + return Category::Normal; + } + + if (grid1 == ECL_GRID_MAINGRID_LGR_NR) { // Main grid + return Category::GlobalToLocal; + } + + return Category::Amalgamated; +} + +bool +Opm::ECLGraph::Impl::NNC:: +isViable(const std::vector& grids, + const int gridID, + const std::size_t cellID) const +{ + using GridIndex = decltype(grids.size()); + const auto gIdx = static_cast(gridID); + + if (gIdx >= grids.size()) { + return false; + } + + const auto& G = grids[gIdx]; + const auto acell = G.activeCell(cellID); + + return (acell >= 0) && (! G.isSubdivided(acell)); +} + +bool +Opm::ECLGraph::Impl::NNC:: +isViable(const std::vector& grids, + const ecl_nnc_type& nnc) const +{ + return (nnc.trans > 0.0) // Omit zero-trans NNCs + && this->isViable(grids, nnc.grid_nr1, nnc.global_index1) + && this->isViable(grids, nnc.grid_nr2, nnc.global_index2); +} + +// ====================================================================== + +Opm::ECLGraph::Impl::Impl(const Path& grid, const Path& init) +{ + const auto G = ECL::loadCase(grid); + auto I = ECL::loadFile(init); + + const auto numGrids = ECL::numGrids(G.get()); + + this->grid_.reserve(numGrids); + this->activeOffset_.reserve(numGrids + 1); + this->activeOffset_.push_back(0); + + for (auto gridID = 0*numGrids; gridID < numGrids; ++gridID) + { + this->grid_.emplace_back(ECL::getGrid(G.get(), gridID), + I.get(), gridID); + + this->activeOffset_.push_back(this->activeOffset_.back() + + this->grid_.back().numCells()); + } + + this->defineNNCs(G.get(), I.get()); +} + +void +Opm::ECLGraph::Impl::assignDataSource(const Path& src) +{ + this->src_ = ECL::loadFile(src); +} + +int +Opm::ECLGraph::Impl:: +numGrids() const +{ + return grid_.size(); +} + +int +Opm::ECLGraph::Impl:: +activeCell(const int gridID, + const std::array& ijk) const +{ + const auto gIdx = + static_castgrid_.size())>(gridID); + + if (gIdx >= this->grid_.size()) { + return -1; + } + + const auto& grid = this->grid_[gIdx]; + + const auto active = grid.activeCell(ijk[0], ijk[1], ijk[2]); + + if ((active < 0) || grid.isSubdivided(active)) { + return -1; + } + + const auto off = static_cast(this->activeOffset_[gIdx]); + + return off + active; +} + +std::size_t +Opm::ECLGraph::Impl::numCells() const +{ + return this->activeOffset_.back(); +} + +std::size_t +Opm::ECLGraph::Impl::numConnections() const +{ + auto nconn = std::size_t{0}; + + for (const auto& G : this->grid_) { + nconn += G.numConnections(); + } + + return nconn + this->nnc_.numConnections(); +} + +std::vector +Opm::ECLGraph::Impl::neighbours() const +{ + auto N = std::vector{}; + + // this->numConnections() includes NNCs. + N.reserve(2 * this->numConnections()); + + { + auto off = this->activeOffset_.begin(); + + for (const auto& G : this->grid_) { + const auto add = static_cast(*off); + + for (const auto& cell : G.neighbours()) { + N.push_back(cell + add); + } + + ++off; + } + } + + { + const auto& nnc = this->nnc_.getNeighbours(); + + N.insert(N.end(), nnc.begin(), nnc.end()); + } + + return N; +} + +std::vector +Opm::ECLGraph::Impl::activePoreVolume() const +{ + auto pvol = std::vector{}; + pvol.reserve(this->numCells()); + + for (const auto& G : this->grid_) { + const auto& pv = G.activePoreVolume(); + + pvol.insert(pvol.end(), pv.begin(), pv.end()); + } + + return pvol; +} + +std::vector +Opm::ECLGraph::Impl:: +flux(const BlackoilPhases::PhaseIndex phase, + const int rptstep) const +{ + if (! ecl_file_has_report_step(this->src_.get(), rptstep)) { + return {}; + } + + const auto vector = this->flowVector(phase); + + auto v = std::vector{}; + + // Recall: this->numConnections() includes NNCs. + const auto totconn = this->numConnections(); + + v.reserve(totconn); + + ecl_file_select_rstblock_report_step(this->src_.get(), rptstep); + + for (const auto& G : this->grid_) { + const auto& q = G.connectionData(this->src_.get(), vector); + + if (q.empty()) { + // Flux vector invalid unless all grids provide this result + // vector data. + return {}; + } + + v.insert(v.end(), q.begin(), q.end()); + } + + if (this->nnc_.numConnections() > 0) { + // Model includes non-neighbouring connections such as faults and/or + // local grid refinement. Extract pertinent flux values for these + // connections. + this->fluxNNC(vector, v); + } + + if (v.size() < totconn) { + // Result vector data not available for NNCs. Entire flux vector is + // invalid. + return {}; + } + + return v; +} + +void +Opm::ECLGraph::Impl::defineNNCs(const ecl_grid_type* G, + const ecl_file_type* init) +{ + for (const auto& nnc : ECL::loadNNC(G, init)) { + this->nnc_.add(this->grid_, this->activeOffset_, nnc); + } +} + +void +Opm::ECLGraph::Impl::fluxNNC(const std::string& vector, + std::vector& flux) const +{ + auto v = std::vector(this->nnc_.numConnections(), 0.0); + auto assigned = std::vector(v.size(), false); + + for (const auto& cat : this->nnc_.allCategories()) { + const auto& rel = this->nnc_.getRelations(cat); + const auto fluxID = vector + rel.fluxID; + + auto gridID = 0; + for (const auto& G : this->grid_) { + const auto& iset = rel.indexSet.getGridCollection(gridID); + + // Must increment grid ID irrespective of early break of + // iteration. + gridID += 1; + + if (iset.empty()) { + // No NNCs for this category in this grid. Skip. + continue; + } + + // Note: Method name is confusing, but does actually do what we + // want here. + const auto q = G.cellData(this->src_.get(), fluxID); + + if (q.empty()) { + // No flux data for this category in this grid. Skip. + continue; + } + + // 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]; + + assigned[ix.neighIdx] = true; + } + } + } + + // NNC flux field is valid if there are no unassigned entries, i.e., if + // there are no 'false' values in the "assigned" record. + const auto valid = assigned.end() == + std::find(assigned.begin(), assigned.end(), false); + + if (valid) { + // This result set (flux) vector is fully supported by at least one + // category of non-Cartesian keywords. Append result to global flux + // vector. + flux.insert(flux.end(), v.begin(), v.end()); + } +} + +std::string +Opm::ECLGraph::Impl:: +flowVector(const BlackoilPhases::PhaseIndex phase) const +{ + const auto vector = std::string("FLR"); // Flow-rate, reservoir + + if (phase == BlackoilPhases::PhaseIndex::Aqua) { + return vector + "WAT"; + } + + if (phase == BlackoilPhases::PhaseIndex::Liquid) { + return vector + "OIL"; + } + + if (phase == BlackoilPhases::PhaseIndex::Vapour) { + return vector + "GAS"; + } + + { + std::ostringstream os; + + os << "Invalid phase index '" << phase << '\''; + + throw std::invalid_argument(os.str()); + } +} + +// ====================================================================== + +Opm::ECLGraph::ECLGraph(ImplPtr pImpl) + : pImpl_(std::move(pImpl)) +{ +} + +Opm::ECLGraph::ECLGraph(ECLGraph&& rhs) + : pImpl_(std::move(rhs.pImpl_)) +{} + +Opm::ECLGraph::~ECLGraph() +{} + +Opm::ECLGraph& +Opm::ECLGraph::operator=(ECLGraph&& rhs) +{ + this->pImpl_ = std::move(rhs.pImpl_); + + return *this; +} + +Opm::ECLGraph +Opm::ECLGraph::load(const Path& grid, const Path& init) +{ + auto pImpl = ImplPtr{new Impl(grid, init)}; + + return { std::move(pImpl) }; +} + +int +Opm::ECLGraph::numGrids() const +{ + return this->pImpl_->numGrids(); +} + +int +Opm::ECLGraph::activeCell(const std::array& ijk, + const int gridID) const +{ + return this->pImpl_->activeCell(gridID, ijk); +} + +void +Opm::ECLGraph::assignFluxDataSource(const Path& src) +{ + this->pImpl_->assignDataSource(src); +} + +std::size_t +Opm::ECLGraph::numCells() const +{ + return this->pImpl_->numCells(); +} + +std::size_t +Opm::ECLGraph::numConnections() const +{ + return this->pImpl_->numConnections(); +} + +std::vector Opm::ECLGraph::neighbours() const +{ + return this->pImpl_->neighbours(); +} + +std::vector Opm::ECLGraph::poreVolume() const +{ + return this->pImpl_->activePoreVolume(); +} + +std::vector +Opm::ECLGraph:: +flux(const BlackoilPhases::PhaseIndex phase, + const int rptstep) const +{ + return this->pImpl_->flux(phase, rptstep); +} diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp new file mode 100644 index 0000000000..9a0570a266 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp @@ -0,0 +1,172 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_ECLGRAPH_HEADER_INCLUDED +#define OPM_ECLGRAPH_HEADER_INCLUDED + +#include + +#include +#include +#include +#include + +#include + +/// \file +/// +/// Facility for extracting active cells and neighbourship relations from +/// on-disk ECLIPSE output, featuring on-demand and cached property loading +/// from backing object (e.g., restart vectors at various time points). + +namespace Opm { + + /// Package an ECLIPSE result set (represented as GRID, INIT, and + /// restart files) in form usable with the computational engine + /// represented by class \code Opm::FlowDiagnostics::ToolBox \endcode. + class ECLGraph + { + public: + using Path = boost::filesystem::path; + + /// Disabled default constructor. + ECLGraph() = delete; + + /// Destructor. + ~ECLGraph(); + + /// Move constructor. + /// + /// \param[in] rhs Graph from which to appropriate resources. + /// Invalid upon return. + ECLGraph(ECLGraph&& rhs); + + /// Disabled copy constructor + ECLGraph(const ECLGraph& rhs) = delete; + + /// Move assignment operator. + /// + /// \param[in] rhs Graph from which to appropriate resources. + /// + /// \return Reference to \code *this \endcode. + ECLGraph& operator=(ECLGraph&& rhs); + + /// Disabled assignment operator. + ECLGraph& operator=(const ECLGraph& rhs) = delete; + + /// Named constructor. + /// + /// \param[in] grid Name or prefix of ECL grid (i.e., .GRID or + /// .EGRID) file. + /// + /// \param[in] init Name of ECL INIT file corresponding to \p grid + /// input. Assumed to provide at least a complete + /// set of pore-volume values (i.e., for all global + /// cells defined in the \p grid). + /// + /// If available in the INIT file, the constructor + /// will also leverage the transmissibility data + /// when constructing the active cell neighbourship + /// table. + /// + /// \return Fully formed ECLIPSE connection graph with property + /// associations. + static ECLGraph + load(const Path& grid, const Path& init); + + /// Assign source object for phase flux calculation. + /// + /// \param[in] src Name of ECL restart file, possibly unified, from + /// which next set of phase fluxes should be retrieved. + void assignFluxDataSource(const Path& src); + + /// Retrieve number of grids. + /// + /// \return The number of LGR grids plus one (the main grid). + int numGrids() const; + + /// Retrieve active cell ID from (I,J,K) tuple in particular grid. + /// + /// \param[in] ijk Cartesian index tuple of particular cell. + /// + /// \param[in] gridID Identity of specific grid to which to relate + /// the (I,J,K) tuple. Use zero (default) for main grid and + /// positive indices for any LGRs. The (I,J,K) indices must be + /// within the ranges implied by the specific grid. + /// + /// \return Active ID (relative to linear, global numbering) of cell + /// (I,J,K) from specified grid. Negative one (-1) if (I,J,K) + /// outside valid range or if the specific cell identified by \p + /// ijk and \p gridID is not actually active. + int activeCell(const std::array& ijk, + const int gridID = 0) const; + + /// Retrieve number of active cells in graph. + std::size_t numCells() const; + + /// Retrive number of connections in graph. + std::size_t numConnections() const; + + /// Retrive 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. + /// + /// 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 poreVolume() const; + + /// Retrive phase flux on all connections defined by \code + /// neighbours() \endcode. + /// + /// \param[in] phase Canonical phase for which to retrive flux. + /// + /// \param[in] rptstep Selected temporal vector. Report-step ID. + /// + /// \return Flux values corresponding to selected phase and report + /// step. Empty if unavailable in the result set (e.g., by querying + /// the gas flux in an oil/water system or if the specified \p + /// occurrence is not reported due to report frequencies or no flux + /// values are output at all). + std::vector + flux(const BlackoilPhases::PhaseIndex phase, + const int rptstep = 0) const; + + private: + /// Implementation class. + class Impl; + + using ImplPtr = std::unique_ptr; + + /// Constructor. Intentially not \c explicit. + ECLGraph(ImplPtr pImpl); + + /// Pointer to implementation. + ImplPtr pImpl_; + }; + +} // namespace Opm + +#endif // OPM_ECLGRAPH_HEADER_INCLUDED diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp new file mode 100644 index 0000000000..81334310b8 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp @@ -0,0 +1,339 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + namespace { + + + + /// RAII class using the ERT block selection stack mechanism + /// to select a report step in a restart file. + struct SelectReportBlock + { + /// \param[in] file ecl file to select block in. + /// \paran[in] report_step sequence number of block to choose + SelectReportBlock(ecl_file_type* file, const int report_step) + : file_(file) + { + if (!ecl_file_has_report_step(file_, report_step)) { + throw std::invalid_argument("Report step " + std::to_string(report_step) + " not found."); + } + ecl_file_push_block(file_); + ecl_file_select_global(file_); + ecl_file_select_rstblock_report_step(file_, report_step); + } + ~SelectReportBlock() + { + ecl_file_pop_block(file_); + } + ecl_file_type* file_; + }; + + + + + /// RAII class using the ERT block selection stack mechanism + /// to select an LGR sub-block. + struct SubSelectLGRBlock + { + /// \param[in] file ecl file to select LGR block in. + /// \paran[in] lgr_index sequence number of block to choose + SubSelectLGRBlock(ecl_file_type* file, const int lgr_index) + : file_(file) + { + ecl_file_push_block(file_); + ecl_file_subselect_block(file_, LGR_KW, lgr_index); + } + ~SubSelectLGRBlock() + { + ecl_file_pop_block(file_); + } + ecl_file_type* file_; + }; + + + + + /// Simple ecl file loading function. + ERT::ert_unique_ptr + load(const boost::filesystem::path& filename) + { + // Read-only, keep open between requests + const auto open_flags = 0; + using FilePtr = ERT::ert_unique_ptr; + FilePtr file(ecl_file_open(filename.generic_string().c_str(), open_flags)); + if (!file) { + std::ostringstream os; + os << "Failed to load ECL File object from '" + << filename.generic_string() << '\''; + throw std::invalid_argument(os.str()); + } + return file; + } + + + + + // --------- Restart file keywords. --------- + + // Note: + // This struct is more complete (containing more fields) + // than currently needed, but we should expect that more + // fields could be needed in the future. + struct INTEHEAD + { + + // Unit codes used in INTEHEAD + enum { + Metric = 1, + Field = 2, + Lab = 3 + }; + + explicit INTEHEAD(const std::vector& v) + : unit (v[INTEHEAD_UNIT_INDEX]) + , nx (v[INTEHEAD_NX_INDEX]) + , ny (v[INTEHEAD_NY_INDEX]) + , nz (v[INTEHEAD_NZ_INDEX]) + , nactive(v[INTEHEAD_NACTIVE_INDEX]) + , iphs (v[INTEHEAD_PHASE_INDEX]) + , nwell (v[INTEHEAD_NWELLS_INDEX]) + , ncwma (v[INTEHEAD_NCWMAX_INDEX]) + , nwgmax (v[INTEHEAD_NWGMAX_INDEX]) + , ngmaxz (v[INTEHEAD_NGMAXZ_INDEX]) + , niwel (v[INTEHEAD_NIWELZ_INDEX]) + , nswel (v[INTEHEAD_NSWELZ_INDEX]) + , nxwel (v[INTEHEAD_NXWELZ_INDEX]) + , nzwel (v[INTEHEAD_NZWELZ_INDEX]) + , nicon (v[INTEHEAD_NICONZ_INDEX]) + , nscon (v[INTEHEAD_NSCONZ_INDEX]) + , nxcon (v[INTEHEAD_NXCONZ_INDEX]) + , nigrpz (v[INTEHEAD_NIGRPZ_INDEX]) + , iday (v[INTEHEAD_DAY_INDEX]) + , imon (v[INTEHEAD_MONTH_INDEX]) + , iyear (v[INTEHEAD_YEAR_INDEX]) + , iprog (v[INTEHEAD_IPROG_INDEX]) + { + } + + int unit; // Unit system. 1:metric, 2:field, 3:lab. + int nx; // Cartesian size i-direction. + int ny; // Cartesian size j-direction. + int nz; // Cartesian size k-direction. + int nactive; // Number of active cells. + int iphs; // Phase. 1:o, 2:w, 3:ow, 4:g, 5:og, 6:wg, 7:owg. + int nwell; // Number of wells. + int ncwma; // Maximum number of completions per well. + int nwgmax; // Maximum number of wells in any group. + int ngmaxz; // Maximum number of groups in field. + int niwel; // Number of elements pr. well in the IWEL array. + int nswel; // Number of elements pr. well in the SWEL array. + int nxwel; // Number of elements pr. well in the XWEL array. + int nzwel; // Number of 8 character words pr. well in ZWEL. + int nicon; // Number of elements pr completion in the ICON array. + int nscon; // Number of elements pr completion in the SCON array. + int nxcon; // Number of elements pr completion in the XCON array. + int nigrpz; // Number of elements pr group in the IGRP array. + int iday; // Report day. + int imon; // Report month. + int iyear; // Report year. + int iprog; // Eclipse program type. 100, 300 or 500. + }; + + + + + // 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 input string with spaces stripped of the right end. + std::string trimSpacesRight(const std::string& s) + { + return std::string(s.begin(), s.begin() + s.find_last_not_of(' ') + 1); + } + + + } // anonymous namespace + + + + + ECLWellSolution::ECLWellSolution(const boost::filesystem::path& restart_filename) + : restart_(load(restart_filename)) + { + } + + + + + + std::vector + ECLWellSolution::solution(const int report_step, + const int num_grids) const + { + SelectReportBlock select(restart_.get(), report_step); + { + // Read well data for global grid. + std::vector all_wd = readWellData(0); + for (int grid_index = 1; grid_index < num_grids; ++grid_index) { + const int lgr_index = grid_index - 1; + SubSelectLGRBlock subselect(restart_.get(), lgr_index); + { + // Read well data for LGR grid. + std::vector wd = readWellData(grid_index); + // Append to set of all well data. + all_wd.insert(all_wd.end(), wd.begin(), wd.end()); + } + } + return all_wd; + } + } + + + + + ecl_kw_type* + ECLWellSolution::getKeyword(const std::string& fieldname) const + { + const int local_occurrence = 0; // This should be correct for all the well-related keywords. + if (ecl_file_get_num_named_kw(restart_.get(), fieldname.c_str()) == 0) { + throw std::runtime_error("Could not find field " + fieldname); + } + return ecl_file_iget_named_kw(restart_.get(), fieldname.c_str(), local_occurrence); + } + + + + + std::vector + ECLWellSolution::loadDoubleField(const std::string& fieldname) const + { + ecl_kw_type* keyword = getKeyword(fieldname); + std::vector field_data; + field_data.resize(ecl_kw_get_size(keyword)); + ecl_kw_get_data_as_double(keyword, field_data.data()); + return field_data; + } + + + + + std::vector + ECLWellSolution::loadIntField(const std::string& fieldname) const + { + ecl_kw_type* keyword = getKeyword(fieldname); + std::vector field_data; + field_data.resize(ecl_kw_get_size(keyword)); + ecl_kw_get_memcpy_int_data(keyword, field_data.data()); + return field_data; + } + + + + + std::vector + ECLWellSolution::loadStringField(const std::string& fieldname) const + { + ecl_kw_type* keyword = getKeyword(fieldname); + std::vector field_data; + const int size = ecl_kw_get_size(keyword); + field_data.resize(size); + for (int pos = 0; pos < size; ++pos) { + field_data[pos] = ecl_kw_iget_char_ptr(keyword, pos); + } + return field_data; + } + + + + + std::vector + ECLWellSolution::readWellData(const int grid_index) const + { + // Note: this function is expected to be called in a context + // where the correct restart block and grid subblock has already + // been selected using the ert block mechanisms. + + // Read header, return if trivial. + INTEHEAD ih(loadIntField(INTEHEAD_KW)); + if (ih.nwell == 0) { + return {}; + } + const double qr_unit = resRateUnit(ih.unit); + + // Read necessary keywords. + auto zwel = loadStringField(ZWEL_KW); + auto iwel = loadIntField(IWEL_KW); + auto icon = loadIntField(ICON_KW); + auto xcon = loadDoubleField("XCON"); + + // Create well data. + std::vector wd(ih.nwell); + for (int well = 0; well < ih.nwell; ++well) { + wd[well].name = trimSpacesRight(zwel[well * ih.nzwel]); + const int ncon = iwel[well * ih.niwel + IWEL_CONNECTIONS_INDEX]; + wd[well].completions.resize(ncon); + for (int comp_index = 0; comp_index < ncon; ++comp_index) { + const int icon_offset = (well*ih.ncwma + comp_index) * ih.nicon; + const int xcon_offset = (well*ih.ncwma + comp_index) * ih.nxcon; + auto& completion = wd[well].completions[comp_index]; + // Note: subtracting 1 from indices (Fortran -> C convention). + completion.grid_index = grid_index; + completion.ijk = { icon[icon_offset + ICON_I_INDEX] - 1, + icon[icon_offset + ICON_J_INDEX] - 1, + icon[icon_offset + ICON_K_INDEX] - 1 }; + // Note: taking the negative input, to get inflow rate. + completion.reservoir_inflow_rate = -unit::convert::from(xcon[xcon_offset + XCON_QR_INDEX], qr_unit); + } + } + return wd; + } + + + + +} // namespace Opm diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp new file mode 100644 index 0000000000..5dfdcf45ce --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp @@ -0,0 +1,79 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_ECLWELLSOLUTION_HEADER_INCLUDED +#define OPM_ECLWELLSOLUTION_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + class ECLWellSolution + { + public: + /// Construct with path to restart file. + explicit ECLWellSolution(const boost::filesystem::path& restart_filename); + + /// Contains the well data extracted from the restart file. + struct WellData + { + std::string name; + struct Completion + { + int grid_index; // 0 for main grid, otherwise LGR grid. + std::array ijk; // Cartesian location in grid. + double reservoir_inflow_rate; // Total fluid rate in SI (m^3/s). + }; + std::vector completions; + }; + + /// Return well solution for given report step. + /// + /// Will throw if required data is not available for the + /// requested step. + std::vector solution(const int report_step, + const int num_grids) const; + + private: + // Types. + using FilePtr = ERT::ert_unique_ptr; + + // Data members. + FilePtr restart_; + + // Methods. + ecl_kw_type* getKeyword(const std::string& fieldname) const; + std::vector loadDoubleField(const std::string& fieldname) const; + std::vector loadIntField(const std::string& fieldname) const; + std::vector loadStringField(const std::string& fieldname) const; + std::vector readWellData(const int grid_index) const; + }; + + +} // namespace Opm + +#endif // OPM_ECLWELLSOLUTION_HEADER_INCLUDED diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opmCore/opm/core/props/BlackoilPhases.hpp b/ThirdParty/custom-opm-flowdiagnostics-applications/opmCore/opm/core/props/BlackoilPhases.hpp new file mode 100644 index 0000000000..a8bff1bf29 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opmCore/opm/core/props/BlackoilPhases.hpp @@ -0,0 +1,78 @@ +/* + Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. + This file is part of the Open Porous Media project (OPM). + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_BLACKOILPHASES_HEADER_INCLUDED +#define OPM_BLACKOILPHASES_HEADER_INCLUDED + + +namespace Opm +{ + + class BlackoilPhases + { + public: + static const int MaxNumPhases = 3; + // enum ComponentIndex { Water = 0, Oil = 1, Gas = 2 }; + enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 }; + + }; + + struct PhaseUsage : public BlackoilPhases + { + int num_phases; + int phase_used[MaxNumPhases]; + int phase_pos[MaxNumPhases]; + }; + + /// Check or assign presence of a formed, free phase. Limited to + /// the 'BlackoilPhases'. + /// + /// Use a std::vector to represent the conditions + /// in an entire model. + class PhasePresence + { + public: + PhasePresence() + : present_(0) + {} + + bool hasFreeWater() const { return present(BlackoilPhases::Aqua ); } + bool hasFreeOil () const { return present(BlackoilPhases::Liquid); } + bool hasFreeGas () const { return present(BlackoilPhases::Vapour); } + + void setFreeWater() { insert(BlackoilPhases::Aqua ); } + void setFreeOil () { insert(BlackoilPhases::Liquid); } + void setFreeGas () { insert(BlackoilPhases::Vapour); } + + bool operator==(const PhasePresence& other) const { return present_ == other.present_; } + bool operator!=(const PhasePresence& other) const { return !this->operator==(other); } + + private: + unsigned char present_; + + bool present(const BlackoilPhases::PhaseIndex i) const + { + return present_ & (1 << i); + } + + void insert(const BlackoilPhases::PhaseIndex i) + { + present_ |= (1 << i); + } + }; + +} // namespace Opm + +#endif // OPM_BLACKOILPHASES_HEADER_INCLUDED \ No newline at end of file diff --git a/ThirdParty/custom-opm-flowdiagnostics-applications/opmCore/opm/core/utility/Units.hpp b/ThirdParty/custom-opm-flowdiagnostics-applications/opmCore/opm/core/utility/Units.hpp new file mode 100644 index 0000000000..33de0ed443 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiagnostics-applications/opmCore/opm/core/utility/Units.hpp @@ -0,0 +1,229 @@ +//=========================================================================== +// +// File: Units.hpp +// +// Created: Thu Jul 2 09:19:08 2009 +// +// Author(s): Halvor M Nilsen +// +// $Date$ +// +// $Revision$ +// +//=========================================================================== + +/* + 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. + */ + +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 + + 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 + /// gives a permeability. We want to set \c k to \f$1\,mD\f$. + /// \code + /// using namespace Opm::unit + /// double k = 0.001*darcy; + /// \endcode + /// We can also use one of the prefixes defined in Opm::prefix + /// \code + /// using namespace Opm::unit + /// 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; } + /// @} + + // -------------------------------------------------------------- + // 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 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 Volume + /// @{ + const double gallon = 231 * cubic(inch); + const double stb = 42 * gallon; + 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; + /// @} + + // -------------------------------------------------------------- + // Standardised constants + // -------------------------------------------------------------- + + /// \name Standardised constant + /// @{ + const double gravity = 9.80665 * meter/square(second); + /// @} + + // -------------------------------------------------------------- + // Derived units and conversions + // -------------------------------------------------------------- + + /// \name Force + /// @{ + const double Newton = kilogram*meter / square(second); // == 1 + 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 Viscosity + /// @{ + const double Pas = Pascal * second; // == 1 + const double Poise = prefix::deci*Pas; + /// @} + + 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; + /// @} + + /** + * 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; + } + + /** + * 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 + + +#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 unit +} // namespace Opm +#endif // OPM_UNITS_HEADER \ No newline at end of file