OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Updated opm-flowdiagnostics-applications to sha: 27adc8964f997d0bb508e254b6ff43ad1e1e76bf

Browse Source
This commit is contained in:
Jacob Støren 2016-12-08 15:19:13 +01:00
parent 77180d91cb
commit 147cb5ebe0
9 changed files with 1537 additions and 478 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/CMakeLists_files.cmake vendored
View File

@ -22,6 +22,7 @@
list (APPEND MAIN_SOURCE_FILES
opm/utility/ECLGraph.cpp
opm/utility/ECLResultData.cpp
opm/utility/ECLWellSolution.cpp
)
@ -34,5 +35,6 @@ list (APPEND EXAMPLE_SOURCE_FILES
list (APPEND PUBLIC_HEADER_FILES
opm/utility/ECLGraph.hpp
opm/utility/ECLResultData.hpp
opm/utility/ECLWellSolution.hpp
)

163
ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/examples/computeToFandTracers.cpp vendored
View File

@ -22,175 +22,16 @@
#include <config.h>
#endif
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include "exampleSetup.hpp"
#include <opm/flowdiagnostics/ConnectivityGraph.hpp>
#include <opm/flowdiagnostics/ConnectionValues.hpp>
#include <opm/flowdiagnostics/Toolbox.hpp>
#include <opm/utility/ECLGraph.hpp>
#include <opm/utility/ECLWellSolution.hpp>
#include <exception>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
#include <boost/filesystem.hpp>
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<std::string>& 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<Opm::ECLWellSolution::WellData>& well_fluxes,
const int step)
{
const auto connGraph = Opm::FlowDiagnostics::
ConnectivityGraph{ static_cast<int>(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=<ecl_case_prefix> step=<report_number>
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<string>("case", "DEFAULT_CASE_NAME");
const path grid = param.has("grid") ? param.get<string>("grid")
: deriveFileName(casename, { ".EGRID", ".FEGRID", ".GRID", ".FGRID" });
const path init = param.has("init") ? param.get<string>("init")
: deriveFileName(casename, { ".INIT", ".FINIT" });
const path restart = param.has("restart") ? param.get<string>("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);
auto fdTool = example::setup(argc, argv);
// Solve for time of flight.
using FDT = Opm::FlowDiagnostics::Toolbox;
std::vector<Opm::FlowDiagnostics::CellSet> start;
auto sol = fdTool.computeInjectionDiagnostics(start);
const auto& tof = sol.fd.timeOfFlight();

227
ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/examples/exampleSetup.hpp vendored Normal file
View File

@ -0,0 +1,227 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_EXAMPLESETUP_HEADER_INCLUDED
#define OPM_EXAMPLESETUP_HEADER_INCLUDED
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/flowdiagnostics/ConnectivityGraph.hpp>
#include <opm/flowdiagnostics/ConnectionValues.hpp>
#include <opm/flowdiagnostics/Toolbox.hpp>
#include <opm/utility/ECLGraph.hpp>
#include <opm/utility/ECLWellSolution.hpp>
#include <exception>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
#include <boost/filesystem.hpp>
namespace example {
inline 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)));
}
inline boost::filesystem::path
deriveFileName(boost::filesystem::path file,
const std::vector<std::string>& 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());
}
inline Opm::FlowDiagnostics::ConnectionValues
extractFluxField(const Opm::ECLGraph& G)
{
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::ECLGraph::PhaseIndex::Aqua ,
Opm::ECLGraph::PhaseIndex::Liquid ,
Opm::ECLGraph::PhaseIndex::Vapour })
{
const auto pflux = G.flux(p);
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;
}
template <class WellFluxes>
Opm::FlowDiagnostics::CellSetValues
extractWellFlows(const Opm::ECLGraph& G,
const WellFluxes& well_fluxes)
{
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);
if (cell_index >= 0) {
inflow.emplace(cell_index, completion.reservoir_inflow_rate);
}
}
}
return inflow;
}
namespace Hack {
inline Opm::FlowDiagnostics::ConnectionValues
convert_flux_to_SI(Opm::FlowDiagnostics::ConnectionValues&& fl)
{
using Co = Opm::FlowDiagnostics::ConnectionValues::ConnID;
using Ph = Opm::FlowDiagnostics::ConnectionValues::PhaseID;
const auto nconn = fl.numConnections();
const auto nphas = fl.numPhases();
for (auto phas = Ph{0}; phas.id < nphas; ++phas.id) {
for (auto conn = Co{0}; conn.id < nconn; ++conn.id) {
fl(conn, phas) /= 86400;
}
}
return fl;
}
}
inline Opm::FlowDiagnostics::Toolbox
initialiseFlowDiagnostics(const Opm::ECLGraph& G)
{
const auto connGraph = Opm::FlowDiagnostics::
ConnectivityGraph{ static_cast<int>(G.numCells()),
G.neighbours() };
// Create the Toolbox.
auto tool = Opm::FlowDiagnostics::Toolbox{ connGraph };
tool.assignPoreVolume(G.poreVolume());
tool.assignConnectionFlux(Hack::convert_flux_to_SI(extractFluxField(G)));
auto wsol = Opm::ECLWellSolution{};
const auto well_fluxes =
wsol.solution(G.rawResultData(), G.numGrids());
tool.assignInflowFlux(extractWellFlows(G, well_fluxes));
return tool;
}
inline auto setup(int argc, char* argv[])
-> decltype(initialiseFlowDiagnostics(std::declval<Opm::ECLGraph>()))
{
// 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<string>("case", "DEFAULT_CASE_NAME");
const path grid = param.has("grid") ? param.get<string>("grid")
: deriveFileName(casename, { ".EGRID", ".FEGRID", ".GRID", ".FGRID" });
const path init = param.has("init") ? param.get<string>("init")
: deriveFileName(casename, { ".INIT", ".FINIT" });
const path restart = param.has("restart") ? param.get<string>("restart")
: deriveFileName(casename, { ".UNRST", ".FUNRST" });
const int step = param.getDefault("step", 0);
// Read graph and fluxes, initialise the toolbox.
auto graph = Opm::ECLGraph::load(grid, init);
graph.assignFluxDataSource(restart);
if (! graph.selectReportStep(step)) {
std::ostringstream os;
os << "Report Step " << step
<< " is Not Available in Result Set '"
<< grid.stem() << '\'';
throw std::domain_error(os.str());
}
return initialiseFlowDiagnostics(graph);
}
} // namespace example
#endif // OPM_EXAMPLESETUP_HEADER_INCLUDED

273
ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp vendored
View File

@ -23,12 +23,15 @@
#endif
#include <opm/utility/ECLGraph.hpp>
#include <opm/utility/ECLResultData.hpp>
#include <algorithm>
#include <array>
#include <cassert>
#include <cstddef>
#include <exception>
#include <initializer_list>
#include <iterator>
#include <map>
#include <memory>
#include <sstream>
@ -37,9 +40,7 @@
#include <boost/filesystem.hpp>
#include <ert/ecl/ecl_file.h>
#include <ert/ecl/ecl_grid.h>
#include <ert/ecl/ecl_kw_magic.h>
#include <ert/ecl/ecl_nnc_export.h>
#include <ert/util/ert_unique_ptr.hpp>
@ -50,7 +51,6 @@
namespace {
namespace ECL {
using GridPtr = ::ERT::ert_unique_ptr<ecl_grid_type, ecl_grid_free>;
using FilePtr = ::ERT::ert_unique_ptr<ecl_file_type, ecl_file_close>;
/// Internalise on-disk representation of ECLIPSE grid.
///
@ -62,13 +62,6 @@ namespace {
/// \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().
@ -110,9 +103,9 @@ namespace {
///
/// \return Vector of pore-volumes for all global cells of \p G.
std::vector<double>
getPVolVector(const ecl_grid_type* G,
const ecl_file_type* init,
const int grid_ID = 0);
getPVolVector(const ecl_grid_type* G,
const ::Opm::ECLResultData& init,
const int grid_ID = 0);
/// Extract non-neighbouring connections from ECLIPSE model
///
@ -124,8 +117,8 @@ namespace {
/// \return Model's non-neighbouring connections, including those
/// between main and local grids.
std::vector<ecl_nnc_type>
loadNNC(const ecl_grid_type* G,
const ecl_file_type* init);
loadNNC(const ecl_grid_type* G,
const ::Opm::ECLResultData& init);
class CartesianGridData
{
@ -140,9 +133,9 @@ namespace {
///
/// \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);
CartesianGridData(const ecl_grid_type* G,
const ::Opm::ECLResultData& init,
const int gridID);
/// Retrieve number of active cells in graph.
std::size_t numCells() const;
@ -192,8 +185,8 @@ namespace {
/// \return Numerical values of result set vector, relative to
/// global cell numbering of this grid.
std::vector<double>
cellData(const ecl_file_type* src,
const std::string& vector) const;
cellData(const ::Opm::ECLResultData& src,
const std::string& vector) const;
/// Retrieve values of result set vector for all Cartesian
/// connections in grid.
@ -206,8 +199,8 @@ namespace {
/// \return Numerical values of result set vector attributed to
/// all of the grid's Cartesian connections.
std::vector<double>
connectionData(const ecl_file_type* src,
const std::string& vector) const;
connectionData(const ::Opm::ECLResultData& src,
const std::string& vector) const;
private:
/// Facility for deriving Cartesian neighbourship in a grid
@ -390,8 +383,8 @@ namespace {
///
/// \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;
bool haveCellData(const ::Opm::ECLResultData& src,
const std::string& keyword) const;
/// Predicate for whether or not a particular result vector is
/// defined on the grid's Cartesian connections.
@ -403,8 +396,8 @@ namespace {
/// \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;
bool haveConnData(const ::Opm::ECLResultData& src,
const std::string& keyword) const;
/// Append directional cell data to global collection of
/// connection data identified by vector name prefix.
@ -419,7 +412,7 @@ namespace {
/// input, collection of values corresponding to any previous
/// directions (preserved), and on output additionally contains
/// the data corresponding to the Cartesian direction \p d.
void connectionData(const ecl_file_type* src,
void connectionData(const ::Opm::ECLResultData& src,
const CartesianCells::Direction d,
const std::string& vector,
std::vector<double>& x) const;
@ -448,7 +441,7 @@ namespace {
///
/// \param[in] init Internalised
void deriveNeighbours(const std::vector<std::size_t>& gcells,
const ecl_file_type* init,
const ::Opm::ECLResultData& init,
const CartesianCells::Direction d);
};
} // namespace ECL
@ -474,9 +467,9 @@ ECL::getGrid(const ecl_grid_type* G, const int gridID)
}
std::vector<double>
ECL::getPVolVector(const ecl_grid_type* G,
const ecl_file_type* init,
const int gridID)
ECL::getPVolVector(const ecl_grid_type* G,
const ::Opm::ECLResultData& init,
const int gridID)
{
auto make_szt = [](const int i)
{
@ -487,14 +480,13 @@ ECL::getPVolVector(const ecl_grid_type* G,
auto pvol = std::vector<double>(nglob, 1.0);
if (ecl_file_has_kw(init, "PORV")) {
auto porv =
ecl_file_iget_named_kw(init, "PORV", gridID);
auto kw = std::string{"PORV"};
assert ((make_szt(ecl_kw_get_size(porv)) == nglob)
&& "Pore-volume must be provided for all global cells");
if (init.haveKeywordData(kw, gridID)) {
pvol = init.keywordData<double>(kw, gridID);
ecl_kw_get_data_as_double(porv, pvol.data());
assert ((pvol.size() == nglob) &&
"Pore-volume must be provided for all global cells");
}
return pvol;
@ -519,28 +511,6 @@ ECL::loadCase(const boost::filesystem::path& grid)
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<std::size_t,3>
ECL::cartesianDimensions(const ecl_grid_type* G)
{
@ -555,8 +525,8 @@ ECL::cartesianDimensions(const ecl_grid_type* G)
}
std::vector<ecl_nnc_type>
ECL::loadNNC(const ecl_grid_type* G,
const ecl_file_type* init)
ECL::loadNNC(const ecl_grid_type* G,
const ::Opm::ECLResultData& init)
{
auto make_szt = [](const int n)
{
@ -570,7 +540,7 @@ ECL::loadNNC(const ecl_grid_type* G,
if (numNNC > 0) {
nncData.resize(numNNC);
ecl_nnc_export(G, init, nncData.data());
ecl_nnc_export(G, init.getRawFilePtr(), nncData.data());
}
return nncData;
@ -783,9 +753,10 @@ CartesianCells::ind2sub(const std::size_t globalCell) const
// ======================================================================
ECL::CartesianGridData::CartesianGridData(const ecl_grid_type* G,
const ecl_file_type* init,
const int gridID)
ECL::CartesianGridData::
CartesianGridData(const ecl_grid_type* G,
const ::Opm::ECLResultData& init,
const int gridID)
: gridID_(gridID)
, cells_ (G, ::ECL::getPVolVector(G, init, gridID_))
{
@ -855,35 +826,31 @@ ECL::CartesianGridData::isSubdivided(const int cellID) const
}
std::vector<double>
ECL::CartesianGridData::cellData(const ecl_file_type* src,
const std::string& vector) const
ECL::CartesianGridData::
cellData(const ::Opm::ECLResultData& 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<double>(ecl_kw_get_size(v));
ecl_kw_get_data_as_double(v, x.data());
auto x = src.keywordData<double>(vector, this->gridID_);
return this->cells_.scatterToGlobal(x);
}
bool
ECL::CartesianGridData::haveCellData(const ecl_file_type* src,
const std::string& vector) const
ECL::CartesianGridData::
haveCellData(const ::Opm::ECLResultData& 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_;
return src.haveKeywordData(vector, this->gridID_);
}
bool
ECL::CartesianGridData::haveConnData(const ecl_file_type* src,
const std::string& vector) const
ECL::CartesianGridData::
haveConnData(const ::Opm::ECLResultData& src,
const std::string& vector) const
{
auto have_data = true;
@ -892,6 +859,7 @@ ECL::CartesianGridData::haveConnData(const ecl_file_type* src,
CartesianCells::Direction::K })
{
const auto vname = this->vectorName(vector, d);
have_data = this->haveCellData(src, vname);
if (! have_data) { break; }
@ -901,8 +869,9 @@ ECL::CartesianGridData::haveConnData(const ecl_file_type* src,
}
std::vector<double>
ECL::CartesianGridData::connectionData(const ecl_file_type* src,
const std::string& vector) const
ECL::CartesianGridData::
connectionData(const ::Opm::ECLResultData& src,
const std::string& vector) const
{
if (! this->haveConnData(src, vector)) {
return {};
@ -922,12 +891,13 @@ ECL::CartesianGridData::connectionData(const ecl_file_type* src,
void
ECL::CartesianGridData::
connectionData(const ecl_file_type* src,
connectionData(const ::Opm::ECLResultData& src,
const CartesianCells::Direction d,
const std::string& vector,
std::vector<double>& x) const
{
const auto v = this->cellData(src, this->vectorName(vector, d));
const auto vname = this->vectorName(vector, d);
const auto v = this->cellData(src, vname);
const auto& cells = this->outCell_.find(d);
@ -955,7 +925,7 @@ vectorName(const std::string& vector,
void
ECL::CartesianGridData::
deriveNeighbours(const std::vector<std::size_t>& gcells,
const ecl_file_type* init,
const ::Opm::ECLResultData& init,
const CartesianCells::Direction d)
{
auto tran = std::string{"TRAN"};
@ -977,7 +947,7 @@ deriveNeighbours(const std::vector<std::size_t>& gcells,
throw std::invalid_argument("Input direction must be (I,J,K)");
}
const auto& T = this->haveCellData(init, tran)
const auto& T = init.haveKeywordData(tran, this->gridID_)
? this->cellData(init, tran)
: std::vector<double>(this->cells_.numGlobalCells(), 1.0);
@ -1079,6 +1049,10 @@ public:
/// strictly positive.
std::vector<double> activePoreVolume() const;
const ::Opm::ECLResultData& rawResultData() const;
bool selectReportStep(const int rptstep) const;
/// Retrive phase flux on all connections defined by \code neighbours()
/// \endcode.
///
@ -1092,8 +1066,7 @@ public:
/// reported due to report frequencies or no flux values are output at
/// all).
std::vector<double>
flux(const BlackoilPhases::PhaseIndex phase,
const int rptstep) const;
flux(const PhaseIndex phase) const;
private:
/// Collection of non-Cartesian neighbourship relations attributed to a
@ -1168,14 +1141,35 @@ private:
/// Map a collection of non-Cartesian neighbourship relations to a
/// specific flux vector identifier.
struct FluxRelation {
class FluxRelation {
public:
explicit FluxRelation(const std::string& fluxID)
: fluxID_(fluxID)
{}
NonNeighKeywordIndexSet& indexSet()
{
return this->indexSet_;
}
const NonNeighKeywordIndexSet& indexSet() const
{
return this->indexSet_;
}
std::string makeKeyword(const std::string& prefix) const
{
return prefix + this->fluxID_;
}
private:
/// Flux vector identifier. Should be one of "N+" for Normal
/// connections, "L+" for GlobalToLocal connections, and "A+"
/// for Amalgamated connections.
std::string fluxID;
std::string fluxID_;
/// Collection of non-Cartesian neighbourship relations.
NonNeighKeywordIndexSet indexSet;
NonNeighKeywordIndexSet indexSet_;
};
/// Constructor.
@ -1267,13 +1261,14 @@ private:
/// Check if connection is viable
///
/// A candidate non-Cartesian connection is viable if both of its
/// endpoints satisfy the viability criterion.
/// A candidate non-Cartesian connection is viable if the associate
/// transmissibility is strictly positive and 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.
/// \return Whether or not the transmissibility is positive and both
/// candidate endpoints satisfy the cell viability criterion.
bool isViable(const std::vector<ECL::CartesianGridData>& grids,
const ecl_nnc_type& nnc) const;
};
@ -1290,7 +1285,7 @@ private:
std::vector<std::size_t> activeOffset_;
/// Current result set.
ECL::FilePtr src_;
std::unique_ptr<ECLResultData> src_;
/// Extract explicit non-neighbouring connections from ECL output.
///
@ -1301,8 +1296,8 @@ private:
/// \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);
void defineNNCs(const ecl_grid_type* G,
const ::Opm::ECLResultData& init);
/// Compute ECL vector basename for particular phase flux.
///
@ -1312,7 +1307,7 @@ private:
/// \return Basename for ECl vector corresponding to particular phase
/// flux.
std::string
flowVector(const BlackoilPhases::PhaseIndex phase) const;
flowVector(const PhaseIndex phase) const;
/// Extract flux values corresponding to particular result set vector
/// for all identified non-neighbouring connections.
@ -1407,7 +1402,10 @@ NNC::add(const std::vector<ECL::CartesianGridData>& grid,
static_cast<std::size_t>(nnc.input_index)
};
this->keywords_[cat].indexSet.add(nnc.grid_nr2, std::move(entry));
auto rel = this->keywords_.find(cat);
if (rel != std::end(this->keywords_)) {
rel->second.indexSet().add(nnc.grid_nr2, std::move(entry));
}
}
std::size_t
@ -1441,13 +1439,13 @@ Opm::ECLGraph::Impl::NNC::makeRelation(const Category cat) const
{
switch (cat) {
case Category::Normal:
return { "N+", {} };
return FluxRelation{ "N+" };
case Category::GlobalToLocal:
return { "L+", {} };
return FluxRelation{ "L+" };
case Category::Amalgamated:
return { "A+", {} };
return FluxRelation{ "A+" };
}
throw std::invalid_argument("Category must be Normal, "
@ -1503,7 +1501,9 @@ isViable(const std::vector<ECL::CartesianGridData>& grids,
Opm::ECLGraph::Impl::Impl(const Path& grid, const Path& init)
{
const auto G = ECL::loadCase(grid);
auto I = ECL::loadFile(init);
auto I = ::Opm::ECLResultData{ init };
I.selectGlobalView();
const auto numGrids = ECL::numGrids(G.get());
@ -1514,19 +1514,19 @@ Opm::ECLGraph::Impl::Impl(const Path& grid, const Path& init)
for (auto gridID = 0*numGrids; gridID < numGrids; ++gridID)
{
this->grid_.emplace_back(ECL::getGrid(G.get(), gridID),
I.get(), gridID);
I, gridID);
this->activeOffset_.push_back(this->activeOffset_.back() +
this->grid_.back().numCells());
}
this->defineNNCs(G.get(), I.get());
this->defineNNCs(G.get(), I);
}
void
Opm::ECLGraph::Impl::assignDataSource(const Path& src)
{
this->src_ = ECL::loadFile(src);
this->src_.reset(new ECLResultData(src));
}
int
@ -1625,15 +1625,21 @@ Opm::ECLGraph::Impl::activePoreVolume() const
return pvol;
}
const ::Opm::ECLResultData&
Opm::ECLGraph::Impl::rawResultData() const
{
return *this->src_;
}
bool Opm::ECLGraph::Impl::selectReportStep(const int rptstep) const
{
return this->src_->selectReportStep(rptstep);
}
std::vector<double>
Opm::ECLGraph::Impl::
flux(const BlackoilPhases::PhaseIndex phase,
const int rptstep) const
flux(const PhaseIndex phase) const
{
if (! ecl_file_has_report_step(this->src_.get(), rptstep)) {
return {};
}
const auto vector = this->flowVector(phase);
auto v = std::vector<double>{};
@ -1643,10 +1649,8 @@ flux(const BlackoilPhases::PhaseIndex phase,
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);
const auto& q = G.connectionData(*this->src_, vector);
if (q.empty()) {
// Flux vector invalid unless all grids provide this result
@ -1674,8 +1678,8 @@ flux(const BlackoilPhases::PhaseIndex phase,
}
void
Opm::ECLGraph::Impl::defineNNCs(const ecl_grid_type* G,
const ecl_file_type* init)
Opm::ECLGraph::Impl::defineNNCs(const ecl_grid_type* G,
const ::Opm::ECLResultData& init)
{
for (const auto& nnc : ECL::loadNNC(G, init)) {
this->nnc_.add(this->grid_, this->activeOffset_, nnc);
@ -1691,11 +1695,11 @@ Opm::ECLGraph::Impl::fluxNNC(const std::string& vector,
for (const auto& cat : this->nnc_.allCategories()) {
const auto& rel = this->nnc_.getRelations(cat);
const auto fluxID = vector + rel.fluxID;
const auto fluxID = rel.makeKeyword(vector);
auto gridID = 0;
for (const auto& G : this->grid_) {
const auto& iset = rel.indexSet.getGridCollection(gridID);
const auto& iset = rel.indexSet().getGridCollection(gridID);
// Must increment grid ID irrespective of early break of
// iteration.
@ -1708,7 +1712,7 @@ Opm::ECLGraph::Impl::fluxNNC(const std::string& vector,
// Note: Method name is confusing, but does actually do what we
// want here.
const auto q = G.cellData(this->src_.get(), fluxID);
const auto q = G.cellData(*this->src_, fluxID);
if (q.empty()) {
// No flux data for this category in this grid. Skip.
@ -1743,26 +1747,27 @@ Opm::ECLGraph::Impl::fluxNNC(const std::string& vector,
std::string
Opm::ECLGraph::Impl::
flowVector(const BlackoilPhases::PhaseIndex phase) const
flowVector(const PhaseIndex phase) const
{
const auto vector = std::string("FLR"); // Flow-rate, reservoir
if (phase == BlackoilPhases::PhaseIndex::Aqua) {
if (phase == PhaseIndex::Aqua) {
return vector + "WAT";
}
if (phase == BlackoilPhases::PhaseIndex::Liquid) {
if (phase == PhaseIndex::Liquid) {
return vector + "OIL";
}
if (phase == BlackoilPhases::PhaseIndex::Vapour) {
if (phase == PhaseIndex::Vapour) {
return vector + "GAS";
}
{
std::ostringstream os;
os << "Invalid phase index '" << phase << '\'';
os << "Invalid phase index '"
<< static_cast<int>(phase) << '\'';
throw std::invalid_argument(os.str());
}
@ -1839,10 +1844,20 @@ std::vector<double> Opm::ECLGraph::poreVolume() const
return this->pImpl_->activePoreVolume();
}
bool Opm::ECLGraph::selectReportStep(const int rptstep) const
{
return this->pImpl_->selectReportStep(rptstep);
}
const Opm::ECLResultData&
Opm::ECLGraph::rawResultData() const
{
return this->pImpl_->rawResultData();
}
std::vector<double>
Opm::ECLGraph::
flux(const BlackoilPhases::PhaseIndex phase,
const int rptstep) const
flux(const PhaseIndex phase) const
{
return this->pImpl_->flux(phase, rptstep);
return this->pImpl_->flux(phase);
}

47
ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp vendored
View File

@ -21,7 +21,7 @@
#ifndef OPM_ECLGRAPH_HEADER_INCLUDED
#define OPM_ECLGRAPH_HEADER_INCLUDED
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/utility/ECLResultData.hpp>
#include <array>
#include <cstddef>
@ -33,8 +33,8 @@
/// \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).
/// on-disk ECLIPSE output, featuring on-demand property loading from
/// backing object (e.g., restart vectors at various time points).
namespace Opm {
@ -121,38 +121,55 @@ namespace Opm {
/// Retrieve number of active cells in graph.
std::size_t numCells() const;
/// Retrive number of connections in graph.
/// Retrieve number of connections in graph.
std::size_t numConnections() const;
/// Retrive neighbourship relations between active cells.
/// Retrieve 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<int> neighbours() const;
/// Retrive static pore-volume values on active cells only.
/// Retrieve static pore-volume values on active cells only.
///
/// Corresponds to the \c PORV vector in the INIT file, possibly
/// restricted to those active cells for which the pore-volume is
/// strictly positive.
std::vector<double> poreVolume() const;
/// Restrict dynamic result set data to single report step.
///
/// This method must be called before calling either flux() or
/// rawResultData().
///
/// \param[in] rptstep Selected temporal vector. Report-step ID.
///
/// \return Whether or not dynamic data for the requested report
/// step exists in the underlying result set identified in method
/// assignFluxDataSource().
bool selectReportStep(const int rptstep) const;
/// Access underlying result set.
///
/// The result set dynamic data corresponds to the most recent call
/// to method selectReportStep().
const ::Opm::ECLResultData& rawResultData() const;
/// Enum for indicating requested phase from the flux() method.
enum class PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
/// Retrive phase flux on all connections defined by \code
/// 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).
/// \return Flux values corresponding to selected phase. Empty if
/// unavailable in the result set (e.g., when querying the gas
/// flux in an oil/water system or if no flux values at all were
/// output to the restart file).
std::vector<double>
flux(const BlackoilPhases::PhaseIndex phase,
const int rptstep = 0) const;
flux(const PhaseIndex phase) const;
private:
/// Implementation class.

883
ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp vendored Normal file
View File

@ -0,0 +1,883 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include <config.h>
#endif
#include <opm/utility/ECLResultData.hpp>
#include <cassert>
#include <ctime>
#include <exception>
#include <initializer_list>
#include <sstream>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
#include <boost/filesystem.hpp>
#include <ert/ecl/ecl_file.h>
#include <ert/ecl/ecl_file_kw.h>
#include <ert/ecl/ecl_file_view.h>
#include <ert/ecl/ecl_grid.h>
#include <ert/ecl/ecl_kw.h>
#include <ert/ecl/ecl_kw_magic.h>
#include <ert/ecl/ecl_nnc_export.h>
#include <ert/ecl/ecl_util.h>
#include <ert/util/ert_unique_ptr.hpp>
/// \file
///
/// Implementation of ECL Result-Set Interface.
namespace {
namespace ECLImpl {
using FilePtr = ::ERT::ert_unique_ptr<ecl_file_type, ecl_file_close>;
namespace Details {
/// Convert vector of elements from one element type to another.
///
/// \tparam Input Element type of input collection.
template <typename Input>
struct Convert
{
/// Convert from one element type to another.
///
/// Implements case of *different* element types.
///
/// \tparam Output Element type of output collection. Must
/// be different from \c Input.
///
/// \param[in] x Input vector.
///
/// \return Result vector (elements of \p x converted to
/// type \p Output).
template <typename Output>
static std::vector<Output>
to(const std::vector<Input>& x, std::false_type)
{
static_assert(! std::is_same<Output, Input>::value,
"Logic Error: Convert<Input>::to<Output>"
" for Output!=Input");
auto result = std::vector<Output>{};
result.reserve(x.size());
for (const auto& xi : x) {
result.emplace_back(xi);
}
return result;
}
/// Convert from one element type to another.
///
/// Implements special case of same element types. This is
/// the identity operator.
///
/// \tparam Output Element type of output collection. Must
/// be the same as \c Input.
///
/// \param[in] x Input vector.
///
/// \return Result vector (copy of \p x).
///
template <typename Output>
static std::vector<Output>
to(const std::vector<Input>& x, std::true_type)
{
static_assert(std::is_same<Output, Input>::value,
"Logic Error: Convert<Input>::to<Output>"
" for Output==Input");
return x;
}
};
/// Retrieve keyword data elements, possibly converted to
/// different destination data type.
///
/// \tparam Output Destination data type. Type of result
/// vector. Assumed to be arithmetic.
///
/// \tparam Input Source data type. Type of keyword data
/// elements. Assumed to be arithmetic.
///
/// \tparam GetData Data element accessor. Callable.
///
/// \param[in] kw Particular ECLIPSE result-set data vector.
///
/// \param[in] extractElements Accessor for keyword data
/// elements. Must support the syntax
///
/// \code
/// extractElements(const ecl_kw_type* kw, Input* x)
/// \endcode
///
/// \return Data elements of keyword as a \code
/// std::vector<Output> \endcode.
template <typename Output, typename Input, class GetData>
std::vector<Output>
getData(const ecl_kw_type* kw, GetData&& extractElements)
{
auto x = std::vector<Input>(ecl_kw_get_size(kw));
extractElements(kw, x.data());
return Convert<Input>::template to<Output>
(x, typename std::is_same<Output, Input>::type());
}
/// Translate ERT type class to keyword element type.
///
/// Primary template.
///
/// \tparam Input Class of ERT keyword data.
template <ecl_type_enum Input>
struct ElementType
{
/// Undefined element type.
using type = void;
};
/// Translate ERT type class to keyword element type.
///
/// Actual element type of \code ECL_INT_TYPE \endcode.
template <>
struct ElementType<ECL_INT_TYPE>
{
/// Element type of ERT integer data.
using type = int;
};
/// Translate ERT type class to keyword element type.
///
/// Actual element type of \code ECL_FLOAT_TYPE \endcode.
template <>
struct ElementType<ECL_FLOAT_TYPE>
{
/// Element type of ERT floating-point (float) data.
using type = float;
};
/// Translate ERT type class to keyword element type.
///
/// Actual element type of \code ECL_DOUBLE_TYPE \endcode.
template <>
struct ElementType<ECL_DOUBLE_TYPE>
{
/// Element type of ERT floating-point (double) data.
using type = double;
};
/// Extract ERT keyword data for various element types.
///
/// Primary template.
///
/// \tparam Input Class of ERT keyword data.
template <ecl_type_enum Input>
struct ExtractKeywordElements;
/// Extract ERT keyword integer data.
template <>
struct ExtractKeywordElements<ECL_INT_TYPE>
{
using EType = ElementType<ECL_INT_TYPE>::type;
/// Function call operator.
///
/// Retrieve actual data elements from ERT keyword of integer
/// (specifically, \c int) type.
///
/// \param[in] kw ERT keyword instance.
///
/// \param[in,out] x Linearised keyword data elements. On
/// input points to memory block of size \code
/// ecl_kw_get_size(kw) * sizeof *x \endcode bytes. On
/// output, those bytes are filled with the actual data
/// values of \p kw.
void operator()(const ecl_kw_type* kw, EType* x) const
{
ecl_kw_get_memcpy_int_data(kw, x);
}
};
/// Extract ERT keyword \c float data.
template <>
struct ExtractKeywordElements<ECL_FLOAT_TYPE>
{
using EType = ElementType<ECL_FLOAT_TYPE>::type;
/// Function call operator.
///
/// Retrieve actual data elements from ERT keyword of
/// floating-point (specifically, \c float) type.
///
/// \param[in] kw ERT keyword instance.
///
/// \param[in,out] x Linearised keyword data elements. On
/// input points to memory block of size \code
/// ecl_kw_get_size(kw) * sizeof *x \endcode bytes. On
/// output, those bytes are filled with the actual data
/// values of \p kw.
void operator()(const ecl_kw_type* kw, EType* x) const
{
ecl_kw_get_memcpy_float_data(kw, x);
}
};
/// Extract ERT keyword \c double data.
template <>
struct ExtractKeywordElements<ECL_DOUBLE_TYPE>
{
using EType = ElementType<ECL_DOUBLE_TYPE>::type;
/// Function call operator.
///
/// Retrieve actual data elements from ERT keyword of
/// floating-point (specifically, \c double) type.
///
/// \param[in] kw ERT keyword instance.
///
/// \param[in,out] x Linearised keyword data elements. On
/// input points to memory block of size \code
/// ecl_kw_get_size(kw) * sizeof *x \endcode bytes. On
/// output, those bytes are filled with the actual data
/// values of \p kw.
void operator()(const ecl_kw_type* kw, EType* x) const
{
ecl_kw_get_memcpy_double_data(kw, x);
}
};
} // namespace Details
/// Extract data elements from ECL keyword.
///
/// Primary template--handles non-string (non-char) types.
///
/// \tparam Input Type of keyword data.
template <ecl_type_enum Input>
struct GetKeywordData
{
/// Retrieve arithmetic (non-string) ECL keyword data as a \code
/// std::vector<T> \endcode for arithmetic types \c T.
///
/// \tparam T Element type of result vector. Must not be a
/// string type.
///
/// \param[in] kw ECL keyword instance. Its element type must
/// be commensurate with \p Input.
///
/// \return Keyword data elements as a \code std::vector<T>
/// \endcode.
template <typename T>
static std::vector<T>
as(const ecl_kw_type* kw, std::false_type)
{
assert (ecl_kw_get_type(kw) == Input);
return Details::getData<
T, typename Details::ElementType<Input>::type
>(kw, Details::ExtractKeywordElements<Input>{});
}
/// Retrieve vector of string data.
///
/// Not supported for non-char keyword element types.
///
/// \tparam T Element type of result vector. Assumed to be
/// \code std::string \endcode.
///
/// \return Empty string vector. String data is unsupported for
/// non-char input types.
template <typename T>
static std::vector<std::string>
as(const ecl_kw_type* /* kw */, std::true_type)
{
return {};
}
};
/// Extract data elements from ECL character-type keyword.
///
/// Only supports accessing elements as a vector of \code
/// std::string \endcode.
template <>
struct GetKeywordData<ECL_CHAR_TYPE>
{
/// Retrieve arithmetic (non-string) ECL keyword data as a \code
/// std::vector<T> \endcode for arithmetic types \c T.
///
/// Not supported for character input types.
///
/// \tparam T Element type of result vector. Must not be a
/// string type.
///
/// \return Empty data vector.
template <typename T>
static std::vector<T>
as(const ecl_kw_type* /* kw */, std::false_type)
{
return {};
}
/// Retrieve vector of string data.
///
/// \tparam T Element type of result vector. Assumed to be
/// \code std::string \endcode.
///
/// \return Keyword data elements as a \code
/// std::vector<std::string> \endcode.
template <typename T>
static std::vector<std::string>
as(const ecl_kw_type* kw, std::true_type)
{
assert (ecl_kw_get_type(kw) == ECL_CHAR_TYPE);
auto result = std::vector<std::string>{};
result.reserve(ecl_kw_get_size(kw));
for (decltype(ecl_kw_get_size(kw))
i = 0, nkw = ecl_kw_get_size(kw);
i < nkw; ++i)
{
result.emplace_back(ecl_kw_iget_char_ptr(kw, i));
}
return result;
}
};
} // namespace ECLImpl
/// Predicate for whether or not a particular path represents a regular
/// file.
///
/// \param[in] p Filesystem element.
///
/// \return Whether or not the element represented by \p p exists and is
/// (or points to, in the case of a symbolic link) a regular file.
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)));
}
/// Derive filesystem element from prefix or filename.
///
/// Pass-through if the input already is a regular file in order to
/// support accessing result-sets other than restart files (e.g., the
/// INIT vectors).
///
/// Fails (throws an exception of type \code std::invalid_argument
/// \endcode) if no valid filesystem element can be derived from the
/// input argument
///
/// \param[in] file Filename or casename prefix.
///
/// \return Filesystem element corresponding to result-set. Either the
/// input \p file itself or, in the case of a casename prefix, the
/// path to a restart file (unified format only).
boost::filesystem::path
deriveResultPath(boost::filesystem::path file)
{
if (isFile(file)) {
return file;
}
for (const auto* ext : { ".UNRST", ".FUNRST" }) {
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());
}
/// Load result-set represented by filesystem element.
///
/// Fails (throws an exception of type \code std::invalid_argument
/// \endcode) if the input filesystem element does not represent a valid
/// result-set resource or if the resource could not be opened (e.g.,
/// due to lack of permission).
///
/// \param[in] rset Fileystem element representing a result-set.
///
/// \return Accessor handle of result-set.
ECLImpl::FilePtr openResultSet(const boost::filesystem::path& rset)
{
// Read-only, keep open between requests
const auto open_flags = 0;
auto F = ECLImpl::FilePtr{
ecl_file_open(rset.generic_string().c_str(), open_flags)
};
if (! F) {
std::ostringstream os;
os << "Failed to load ECL Result object from '"
<< rset.generic_string() << '\'';
throw std::invalid_argument(os.str());
}
return F;
}
/// Retrieve first keyword from file
///
/// Useful in order to identify a result-set as either a unified restart
/// file or some other result-set type.
///
/// \param[in] file Raw result-set.
///
/// \return First keyword in result-set represented by \p file.
std::string firstFileKeyword(const ecl_file_type* file)
{
// Note: ecl_file_get_global_view() does not modify its input.
auto* globView =
ecl_file_get_global_view(const_cast<ecl_file_type*>(file));
return ecl_kw_get_header(ecl_file_view_iget_kw(globView, 0));
}
} // namespace Anonymous
/// Engine powering implementation of \c ECLResultData interface
class Opm::ECLResultData::Impl
{
public:
/// Constructor
///
/// \param[in] rset Filesystem element or casename prefix representing
/// an ECL result-set.
Impl(Path rset);
/// Copy constructor.
///
/// \param[in] rhs Object from which to construct new \c Impl instance.
Impl(const Impl& rhs);
/// Move constructor.
///
/// \param[in,out] rhs Object from which to constructo new \c Impl
/// instance. Underlying result-set accessor is null upon return.
Impl(Impl&& rhs);
/// Access the underlying ERT representation of the result-set.
///
/// This is a hole in the interface that exists to be able to access
/// ERT's internal data structures for non-neighbouring connections
/// (i.e., faults and/or connections between main grid and LGRs). See
/// function ecl_nnc_export() in the ERT.
///
/// Handle with care.
///
/// \return Handle to underlying ERT representation of result-set.
const ecl_file_type* getRawFilePtr() const;
/// Reset the object's internal view of the result-set to encompass the
/// entirety of the underlying file's result vectors.
///
/// This is mostly useful in the case of accessing static result-sets
/// such as INIT files.
///
/// \return Whether or not generating the global view succeeded.
bool selectGlobalView();
/// Select a result-set view that corresponds to a single report step.
///
/// This is needed when working with dynamic restart data.
///
/// \param[in] step Report step number.
///
/// \return Whether or not selecting the report step succeeded. The
/// typical failure case is the report step not being available
/// in the result-set.
bool selectReportStep(const int step);
/// Query current result-set view for availability of particular named
/// result vector in particular enumerated grid.
///
/// \param[in] vector Named result vector for which to query data
/// availability.
///
/// \param[in] gridID Identity of specific grid for which to query data
/// availability.
///
/// \return Whether or not keyword data for the named result vector is
/// available in the specific grid.
bool haveKeywordData(const std::string& vector,
const int gridID) const;
/// Retrieve current result-set view's data values for particular named
/// result vector in particular enumerated grid.
///
/// Will fail (throw an exception of type std::invalid_argument) unless
/// the requested keyword data is available in the specific grid in the
/// current active view.
///
/// \tparam T Element type of return value. The underlying keyword data
/// will be converted to this type if needed and possible. Note that
/// some type combinations do not make sense. It is, for instance,
/// not possible to retrieve keyword values of an underlying
/// arithmetic type in the form of a \code std::vector<std::string>
/// \endcode. Similarly, we cannot access underlying character data
/// through elements of an arithmetic type (e.g., \code
/// std::vector<double> \endcode.)
///
/// \param[in] vector Named result vector for which to retrieve
/// keyword data.
///
/// \param[in] gridID Identity of specific grid for which to
/// retrieve keyword data.
///
/// \return Keyword data values. Empty if type conversion fails.
template <typename T>
std::vector<T>
keywordData(const std::string& vector,
const int gridID) const;
private:
/// RAII class to select a sub-block pertaining to a particular grid
/// within the current active result-set view.
///
/// An object of this type restricts the view to a particular grid in
/// the constructor and restores the original active view in the
/// destructor.
class Restrict
{
public:
/// Constructor.
///
/// Restricts active result-set view to particular grid.
///
/// \param[in] host Result-set host object that maintains the
/// current active view.
///
/// \param[in] gridID Identity of specific grid to which to restrict
/// the current active result-set view.
Restrict(const Impl& host, const int gridID)
: host_(host)
, save_(host_.activeBlock_)
{
if (gridID == ECL_GRID_MAINGRID_LGR_NR) {
const auto& start = this->host_.mainGridStart();
this->host_.activeBlock_ =
ecl_file_view_add_blockview2(this->host_.activeBlock_,
start.c_str(), LGR_KW, 0);
}
else if (gridID > ECL_GRID_MAINGRID_LGR_NR) {
this->host_.activeBlock_ =
ecl_file_view_add_blockview2(this->host_.activeBlock_,
LGR_KW, LGR_KW, gridID - 1);
}
}
/// Destructor.
///
/// Restores original active result-set view (i.e., widens view to
/// encompass all grids).
~Restrict()
{
this->host_.activeBlock_ = this->save_;
}
private:
/// Host object that maintains the active result-set view.
const Impl& host_;
/// Saved original active view from host (prior to restricting view
/// to single grid ID).
ecl_file_view_type* save_;
};
/// Casename prefix. Mostly to implement copy ctor.
const Path prefix_;
/// Active result-set.
ECLImpl::FilePtr result_;
/// First keyword in result-set (\c result_). Needed to identify start
/// of main grid's section within a view.
std::string firstKeyword_;
/// Current active result-set view.
mutable ecl_file_view_type* activeBlock_{ nullptr };
/// Support for passing \code *this \endcode to ERT functions that
/// require an \c ecl_file_type, particularly the function that selects
/// the file's global view--ecl_file_get_global_view().
///
/// Note mutable return type. This is okay because the relevant
/// functions don't modify their inputs.
operator ecl_file_type*() const;
/// Support for passing \code *this \endcode to ERT functions that
/// require an \c ecl_file_view_type.
operator const ecl_file_view_type*() const;
/// Retrieve result-set keyword that identifies beginning of main grid's
/// result vectors.
const std::string& mainGridStart() const;
};
Opm::ECLResultData::Impl::Impl(Path prefix)
: prefix_ (std::move(prefix))
, result_ (openResultSet(deriveResultPath(prefix_)))
, firstKeyword_(firstFileKeyword(result_.get()))
{}
Opm::ECLResultData::Impl::Impl(const Impl& rhs)
: prefix_ (rhs.prefix_)
, result_ (openResultSet(deriveResultPath(prefix_)))
, firstKeyword_(firstFileKeyword(result_.get()))
{}
Opm::ECLResultData::Impl::Impl(Impl&& rhs)
: prefix_ (std::move(rhs.prefix_))
, result_ (std::move(rhs.result_))
, firstKeyword_(std::move(rhs.firstKeyword_))
{}
const ecl_file_type*
Opm::ECLResultData::Impl::getRawFilePtr() const
{
return this->result_.get();
}
bool Opm::ECLResultData::Impl::selectGlobalView()
{
this->activeBlock_ = ecl_file_get_global_view(*this);
return this->activeBlock_ != nullptr;
}
bool Opm::ECLResultData::Impl::selectReportStep(const int step)
{
if (! ecl_file_has_report_step(*this, step)) {
return false;
}
if (auto* globView = ecl_file_get_global_view(*this)) {
// Ignore sequence numbers, dates, and simulation time.
const auto seqnum = -1;
const auto dates = static_cast<std::size_t>(-1);
const auto simdays = -1.0;
this->activeBlock_ =
ecl_file_view_add_restart_view(globView, seqnum,
step, dates, simdays);
return this->activeBlock_ != nullptr;
}
return false;
}
bool Opm::ECLResultData::Impl::
haveKeywordData(const std::string& vector, const int gridID) const
{
assert ((gridID >= 0) && "Grid IDs must be non-negative");
// Note: Non-trivial dtor. Compiler can't ignore object.
const auto block = Restrict{ *this, gridID };
const auto count =
ecl_file_view_get_num_named_kw(*this, vector.c_str());
return count > 0;
}
namespace Opm {
template <typename T>
std::vector<T>
ECLResultData::Impl::keywordData(const std::string& vector,
const int gridID) const
{
if (! this->haveKeywordData(vector, gridID)) {
std::ostringstream os;
os << "Cannot Access Non-Existent Keyword Data Pair ("
<< vector << ", " << gridID << ')';
throw std::invalid_argument(os.str());
}
// Note: Non-trivial dtor. Compiler can't ignore object.
const auto block = Restrict{ *this, gridID };
const auto occurrence = 0;
const auto* kw =
ecl_file_view_iget_named_kw(*this, vector.c_str(),
occurrence);
assert ((kw != nullptr) &&
"Logic Error In Data Availability Check");
// Whether or not caller requests a vector<string>.
const auto makeStringVector =
typename std::is_same<T, std::string>::type{};
switch (ecl_kw_get_type(kw)) {
case ECL_CHAR_TYPE:
return ECLImpl::GetKeywordData<ECL_CHAR_TYPE>::
as<T>(kw, makeStringVector);
case ECL_INT_TYPE:
return ECLImpl::GetKeywordData<ECL_INT_TYPE>::
as<T>(kw, makeStringVector);
case ECL_FLOAT_TYPE:
return ECLImpl::GetKeywordData<ECL_FLOAT_TYPE>::
as<T>(kw, makeStringVector);
case ECL_DOUBLE_TYPE:
return ECLImpl::GetKeywordData<ECL_DOUBLE_TYPE>::
as<T>(kw, makeStringVector);
default:
// No operator exists for this type. Return empty.
return {};
}
}
} // namespace Opm
Opm::ECLResultData::Impl::operator ecl_file_type*() const
{
return this->result_.get();
}
Opm::ECLResultData::Impl::operator const ecl_file_view_type*() const
{
return this->activeBlock_;
}
const std::string&
Opm::ECLResultData::Impl::mainGridStart() const
{
return this->firstKeyword_;
}
// ======================================================================
// Implementation of class Opm::ECLResultData Below Separator
// ======================================================================
Opm::ECLResultData::ECLResultData(const Path& prefix)
: pImpl_(new Impl(prefix))
{}
Opm::ECLResultData::ECLResultData(const ECLResultData& rhs)
: pImpl_(new Impl(*rhs.pImpl_))
{}
Opm::ECLResultData::ECLResultData(ECLResultData&& rhs)
: pImpl_(std::move(rhs.pImpl_))
{}
Opm::ECLResultData&
Opm::ECLResultData::operator=(const ECLResultData& rhs)
{
this->pImpl_.reset(new Impl(*rhs.pImpl_));
return *this;
}
Opm::ECLResultData&
Opm::ECLResultData::operator=(ECLResultData&& rhs)
{
this->pImpl_ = std::move(rhs.pImpl_);
return *this;
}
Opm::ECLResultData::~ECLResultData()
{}
const ecl_file_type*
Opm::ECLResultData::getRawFilePtr() const
{
return this->pImpl_->getRawFilePtr();
}
bool Opm::ECLResultData::selectGlobalView()
{
return this->pImpl_->selectGlobalView();
}
bool Opm::ECLResultData::selectReportStep(const int step)
{
return this->pImpl_->selectReportStep(step);
}
bool
Opm::ECLResultData::
haveKeywordData(const std::string& vector, const int gridID) const
{
return this->pImpl_->haveKeywordData(vector, gridID);
}
namespace Opm {
template <typename T>
std::vector<T>
ECLResultData::keywordData(const std::string& vector,
const int gridID) const
{
return this->pImpl_->template keywordData<T>(vector, gridID);
}
// Explicit instantiations for those types we care about.
template std::vector<std::string>
ECLResultData::keywordData<std::string>(const std::string& vector,
const int gridID) const;
template std::vector<int>
ECLResultData::keywordData<int>(const std::string& vector,
const int gridID) const;
template std::vector<double>
ECLResultData::keywordData<double>(const std::string& vector,
const int gridID) const;
} // namespace Opm::ECL

178
ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLResultData.hpp vendored Normal file
View File

@ -0,0 +1,178 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ECLRESTARTDATA_HEADER_INCLUDED
#define OPM_ECLRESTARTDATA_HEADER_INCLUDED
#include <memory>
#include <string>
#include <vector>
#include <boost/filesystem/path.hpp>
/// \file
///
/// Interface to direct result-set data vector operations.
/// Forward-declaration of ERT's representation of an ECLIPSE result file.
///
/// This is a hole in the insulation between the interface and the
/// underlying implementation of class ECLResultData.
extern "C" {
typedef struct ecl_file_struct ecl_file_type;
} // extern "C"
namespace Opm {
/// Representation of an ECLIPSE result-set.
///
/// This class is aware of the internal structure of ECLIPSE restart
/// files and may restrict its operation to a single report step. The
/// class furthermore knows about sub-blocks corresponding to main or
/// local grids within a report step and queries only those objects that
/// pertain to a single grid at a time.
///
/// Note: The client must select a view of the result-set before
/// accessing any vectors within the set.
class ECLResultData
{
public:
using Path = boost::filesystem::path;
/// Default constructor disabled.
ECLResultData() = delete;
/// Constructor.
///
/// \param[in] casePrefix Name or prefix of ECL result data.
explicit ECLResultData(const Path& casePrefix);
/// Copy constructor.
///
/// \param[in] rhs Object from which to construct new instance.
ECLResultData(const ECLResultData& rhs);
/// Move constructor.
///
/// \param[in,out] rhs Object from which to construct new instance.
/// Its internal implementation will be subsumed into the new
/// object.
ECLResultData(ECLResultData&& rhs);
/// Assignment operator.
///
/// \param[in] rhs Object from which to assign new values to current
/// instance.
///
/// \return \code *this \endcode.
ECLResultData& operator=(const ECLResultData& rhs);
/// Move assignment operator.
///
/// \param[in,out] Object from which to assign new instance values.
/// Its internal implementation will be subsumed into this
/// instance.
///
/// \return \code *this \endcode.
ECLResultData& operator=(ECLResultData&& rhs);
/// Destructor.
~ECLResultData();
/// Access the underlying ERT representation of the result-set.
///
/// This is essentially a hole in the interface that is intended to
/// support a few very specialised uses. Handle with care.
///
/// \return Handle to underlying ERT representation of result-set.
const ecl_file_type* getRawFilePtr() const;
/// Reset the object's internal view of the result-set to encompass
/// the entirety of the underlying file's result vectors.
///
/// This is mostly useful in the case of accessing static result
/// sets such as INIT files.
///
/// \return Whether or not generating the global view succeeded.
bool selectGlobalView();
/// Select a result-set view that corresponds to a single report
/// step.
///
/// This is needed when working with dynamic restart data.
///
/// \param[in] step Report step number.
///
/// \return Whether or not selecting the report step succeeded. The
/// typical failure case is the report step not being available
/// in the result-set.
bool selectReportStep(const int step);
/// Query current result-set view for availability of particular
/// named result vector in particular enumerated grid.
///
/// \param[in] vector Named result vector for which to query data
/// availability.
///
/// \param[in] gridID Identity of specific grid for which to query
/// data availability.
///
/// \return Whether or not keyword data for the named result vector
/// is available in the specific grid.
bool haveKeywordData(const std::string& vector,
const int gridID) const;
/// Retrieve current result-set view's data values for particular
/// named result vector in particular enumerated grid.
///
/// Will fail (throw an exception of type std::invalid_argument)
/// unless the requested keyword data is available in the specific
/// grid in the current active view.
///
/// \tparam T Element type of return value. The underlying keyword
/// data will be converted to this type if needed and possible.
/// Note that some type combinations do not make sense. It is,
/// for instance, not possible to retrieve keyword values of an
/// underlying arithmetic type in the form of a \code
/// std::vector<std::string> \endcode. Similarly, we cannot
/// access underlying character data through elements of an
/// arithmetic type (e.g., \code std::vector<double> \endcode.)
///
/// \param[in] vector Named result vector for which to retrieve
/// keyword data.
///
/// \param[in] gridID Identity of specific grid for which to
/// retrieve keyword data.
///
/// \return Keyword data values. Empty if type conversion fails.
template <typename T>
std::vector<T>
keywordData(const std::string& vector,
const int gridID) const;
private:
class Impl;
std::unique_ptr<Impl> pImpl_;
};
} // namespace Opm
#endif // OPM_ECLRESTARTDATA_HEADER_INCLUDED

212
ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp vendored
View File

@ -23,9 +23,11 @@
#endif
#include <opm/utility/ECLWellSolution.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/utility/ECLResultData.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <ert/ecl/ecl_kw_magic.h>
#include <ert/ecl_well/well_const.h>
#include <cmath>
#include <stdexcept>
#include <sstream>
@ -34,76 +36,6 @@ 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<ecl_file_type, ecl_file_close>
load(const boost::filesystem::path& filename)
{
// Read-only, keep open between requests
const auto open_flags = 0;
using FilePtr = ERT::ert_unique_ptr<ecl_file_type, ecl_file_close>;
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:
@ -197,13 +129,22 @@ namespace Opm
}
// Constants not provided by ert.
enum { XWEL_RESV_INDEX = 4 };
enum { IWEL_TYPE_PRODUCER = 1 };
} // anonymous namespace
ECLWellSolution::ECLWellSolution(const boost::filesystem::path& restart_filename)
: restart_(load(restart_filename))
ECLWellSolution::ECLWellSolution(const double rate_threshold,
const bool disallow_crossflow)
: rate_threshold_(rate_threshold)
, disallow_crossflow_(disallow_crossflow)
{
}
@ -212,115 +153,62 @@ namespace Opm
std::vector<ECLWellSolution::WellData>
ECLWellSolution::solution(const int report_step,
ECLWellSolution::solution(const ECLResultData& restart,
const int num_grids) const
{
SelectReportBlock select(restart_.get(), report_step);
{
// Read well data for global grid.
std::vector<WellData> 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<WellData> wd = readWellData(grid_index);
// Append to set of all well data.
all_wd.insert(all_wd.end(), wd.begin(), wd.end());
}
}
return all_wd;
// Read well data for global grid.
std::vector<WellData> all_wd{};
for (int grid_index = 0; grid_index < num_grids; ++grid_index) {
std::vector<WellData> wd = readWellData(restart, grid_index);
// Append to set of all well data.
all_wd.insert(all_wd.end(), wd.begin(), wd.end());
}
}
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<double>
ECLWellSolution::loadDoubleField(const std::string& fieldname) const
{
ecl_kw_type* keyword = getKeyword(fieldname);
std::vector<double> 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<int>
ECLWellSolution::loadIntField(const std::string& fieldname) const
{
ecl_kw_type* keyword = getKeyword(fieldname);
std::vector<int> 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<std::string>
ECLWellSolution::loadStringField(const std::string& fieldname) const
{
ecl_kw_type* keyword = getKeyword(fieldname);
std::vector<std::string> 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;
return all_wd;
}
std::vector<ECLWellSolution::WellData>
ECLWellSolution::readWellData(const int grid_index) const
ECLWellSolution::readWellData(const ECLResultData& restart, 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.
// where the correct restart block using the ert block mechanisms.
// Read header, return if trivial.
INTEHEAD ih(loadIntField(INTEHEAD_KW));
INTEHEAD ih(restart.keywordData<int>(INTEHEAD_KW, grid_index));
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");
auto zwel = restart.keywordData<std::string>(ZWEL_KW, grid_index);
auto iwel = restart.keywordData<int> (IWEL_KW, grid_index);
auto xwel = restart.keywordData<double> ("XWEL" , grid_index);
auto icon = restart.keywordData<int> (ICON_KW, grid_index);
auto xcon = restart.keywordData<double> ("XCON" , grid_index);
// Create well data.
std::vector<WellData> wd(ih.nwell);
std::vector<WellData> wd_vec;
wd_vec.reserve(ih.nwell);
for (int well = 0; well < ih.nwell; ++well) {
wd[well].name = trimSpacesRight(zwel[well * ih.nzwel]);
// Skip if total rate below threshold (for wells that are
// shut or stopped for example).
const double well_reservoir_inflow_rate = -unit::convert::from(xwel[well * ih.nxwel + XWEL_RESV_INDEX], qr_unit);
if (std::fabs(well_reservoir_inflow_rate) < rate_threshold_) {
continue;
}
// Otherwise: add data for this well.
WellData wd;
wd.name = trimSpacesRight(zwel[well * ih.nzwel]);
const int ncon = iwel[well * ih.niwel + IWEL_CONNECTIONS_INDEX];
wd[well].completions.resize(ncon);
const bool is_producer = (iwel[well * ih.niwel + IWEL_TYPE_INDEX] == IWEL_TYPE_PRODUCER);
wd.completions.reserve(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];
WellData::Completion completion;
// Note: subtracting 1 from indices (Fortran -> C convention).
completion.grid_index = grid_index;
completion.ijk = { icon[icon_offset + ICON_I_INDEX] - 1,
@ -328,9 +216,19 @@ namespace Opm
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);
if (disallow_crossflow_) {
// Add completion only if not cross-flowing (injecting producer or producing injector).
if (completion.reservoir_inflow_rate < 0.0 == is_producer) {
wd.completions.push_back(completion);
}
} else {
// Always add completion.
wd.completions.push_back(completion);
}
}
wd_vec.push_back(wd);
}
return wd;
return wd_vec;
}

30
ThirdParty/custom-opm-flowdiagnostics-applications/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp vendored
View File

@ -21,9 +21,6 @@
#ifndef OPM_ECLWELLSOLUTION_HEADER_INCLUDED
#define OPM_ECLWELLSOLUTION_HEADER_INCLUDED
#include <ert/ecl/ecl_file.h>
#include <ert/util/ert_unique_ptr.hpp>
#include <boost/filesystem.hpp>
#include <array>
#include <string>
#include <utility>
@ -32,16 +29,22 @@
namespace Opm
{
class ECLResultData;
class ECLWellSolution
{
public:
/// Construct with path to restart file.
explicit ECLWellSolution(const boost::filesystem::path& restart_filename);
/// Constructor.
/// \param[in] rate_threshold a well will be ignored if its total RESV rate is less than this (m^3/s)
/// \param[in] disallow_crossflow if true, injecting perforations of production wells and vice versa will be ignored
explicit ECLWellSolution(const double rate_threshold = 1e-14,
const bool disallow_crossflow = true);
/// Contains the well data extracted from the restart file.
struct WellData
{
std::string name;
bool is_injector_well;
struct Completion
{
int grid_index; // 0 for main grid, otherwise LGR grid.
@ -51,26 +54,21 @@ namespace Opm
std::vector<Completion> completions;
};
/// Return well solution for given report step.
/// Return well solution for pre-selected report step
///
/// Will throw if required data is not available for the
/// requested step.
std::vector<WellData> solution(const int report_step,
std::vector<WellData> solution(const ECLResultData& restart,
const int num_grids) const;
private:
// Types.
using FilePtr = ERT::ert_unique_ptr<ecl_file_type, ecl_file_close>;
// Data members.
FilePtr restart_;
double rate_threshold_;
bool disallow_crossflow_;
// Methods.
ecl_kw_type* getKeyword(const std::string& fieldname) const;
std::vector<double> loadDoubleField(const std::string& fieldname) const;
std::vector<int> loadIntField(const std::string& fieldname) const;
std::vector<std::string> loadStringField(const std::string& fieldname) const;
std::vector<WellData> readWellData(const int grid_index) const;
std::vector<WellData> readWellData(const ECLResultData& restart,
const int grid_index) const;
};