From 0d35d64aea498181da394fcb05b03be3a6d1c1a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 22 Jun 2021 17:12:45 +0200 Subject: [PATCH] Load Analytic Aquifer Objects from Restart File This commit adds support for creating the analytic aquifer objects Opm::AquiferCT Opm::Aquifetp Opm::Aquancon Opm::AquiferConfig from information stored in the restart vectors {I,S,X}AAQ {I,S}CAQ We add a new helper class Opm::RestartIO::RstAquifer which contain the same data members as the '*_data' structures of the analytic aquifer objects. Those analytic aquifer objects then get friendship from the '*_data' structures in order to assign the private members from the corresponding restart information. We finally add a gateway to EclipseState that consumes an RstAquifer instance and overwrites the internal AquiferConfig object when the restarted run contains analytic aquifers. Update RstState constructor API to meet requirements of RstAquifer, notably by adding an EclipseGrid parameter. That in turn is needed by the RstAquifer to translate connection (I,J,K) tuples to active cell IDs. Note that if an analytic aquifer does not have any connections then this facility will not load said aquifer. That may be a problem when plotting summary curves, but we will address the issue later if it comes up. --- CMakeLists_files.cmake | 2 + opm/io/eclipse/RestartFileView.hpp | 2 + opm/io/eclipse/rst/aquifer.hpp | 132 ++++ opm/io/eclipse/rst/state.hpp | 16 +- .../eclipse/EclipseState/Aquifer/Aquancon.hpp | 7 +- .../EclipseState/Aquifer/AquiferCT.hpp | 9 + .../EclipseState/Aquifer/AquiferConfig.hpp | 16 +- .../eclipse/EclipseState/Aquifer/Aquifetp.hpp | 9 + .../eclipse/EclipseState/EclipseState.hpp | 13 +- src/opm/io/eclipse/RestartFileView.cpp | 37 +- src/opm/io/eclipse/rst/aquifer.cpp | 683 ++++++++++++++++++ src/opm/io/eclipse/rst/state.cpp | 26 +- .../eclipse/EclipseState/Aquifer/Aquancon.cpp | 43 +- .../EclipseState/Aquifer/AquiferCT.cpp | 51 +- .../EclipseState/Aquifer/AquiferConfig.cpp | 8 + .../eclipse/EclipseState/Aquifer/Aquifetp.cpp | 41 ++ .../eclipse/EclipseState/EclipseState.cpp | 8 +- 17 files changed, 1070 insertions(+), 33 deletions(-) create mode 100644 opm/io/eclipse/rst/aquifer.hpp create mode 100644 src/opm/io/eclipse/rst/aquifer.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 4d427bb8c..c61784a57 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -284,6 +284,7 @@ if(ENABLE_ECL_OUTPUT) src/opm/io/eclipse/ExtSmryOutput.cpp src/opm/io/eclipse/RestartFileView.cpp src/opm/io/eclipse/SummaryNode.cpp + src/opm/io/eclipse/rst/aquifer.cpp src/opm/io/eclipse/rst/connection.cpp src/opm/io/eclipse/rst/group.cpp src/opm/io/eclipse/rst/header.cpp @@ -864,6 +865,7 @@ if(ENABLE_ECL_OUTPUT) opm/io/eclipse/ExtSmryOutput.hpp opm/io/eclipse/RestartFileView.hpp opm/io/eclipse/SummaryNode.hpp + opm/io/eclipse/rst/aquifer.hpp opm/io/eclipse/rst/connection.hpp opm/io/eclipse/rst/group.hpp opm/io/eclipse/rst/header.hpp diff --git a/opm/io/eclipse/RestartFileView.hpp b/opm/io/eclipse/RestartFileView.hpp index 44ad59a17..7d39d27ce 100644 --- a/opm/io/eclipse/RestartFileView.hpp +++ b/opm/io/eclipse/RestartFileView.hpp @@ -57,6 +57,8 @@ public: getKeyword(const std::string& vector, const int occurrence = 0) const; const std::vector& intehead() const; + const std::vector& logihead() const; + const std::vector& doubhead() const; private: class Implementation; diff --git a/opm/io/eclipse/rst/aquifer.hpp b/opm/io/eclipse/rst/aquifer.hpp new file mode 100644 index 000000000..13be5b83d --- /dev/null +++ b/opm/io/eclipse/rst/aquifer.hpp @@ -0,0 +1,132 @@ +/* + Copyright 2021 Equinor 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_RESTART_AQUIFER_HPP +#define OPM_RESTART_AQUIFER_HPP + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { + class AquiferConfig; + class EclipseGrid; + class UnitSystem; +} // Opm + +namespace Opm { namespace EclIO { + class RestartFileView; +}} // Opm::EclIO + +namespace Opm { namespace RestartIO { + + class RstAquifer + { + public: + struct CarterTracy { + int aquiferID{}; + int inftableID{}; + int pvttableID{}; + + double porosity{}; + double datum_depth{}; + double total_compr{}; + double inner_radius{}; + double permeability{}; + double thickness{}; + double angle_fraction{}; + double initial_pressure{}; + + double time_constant{}; + double influx_constant{}; + double water_density{}; + double water_viscosity{}; + }; + + struct Fetkovich { + int aquiferID{}; + int pvttableID{}; + + double prod_index{}; + double total_compr{}; + double initial_watvolume{}; + double datum_depth{}; + + double initial_pressure{}; + double time_constant{}; + }; + + class Connections { + public: + struct Cell { + std::size_t global_index; + double influx_coeff; + double effective_facearea; + FaceDir::DirEnum face_dir; + }; + + const std::vector& cells() const + { + return this->cells_; + } + + void reserve(const std::vector::size_type cpty) + { + this->cells_.reserve(cpty); + } + + template + void emplace_back(Args&&... args) + { + this->cells_.push_back(Cell { std::forward(args)... }); + } + + private: + std::vector cells_{}; + }; + + explicit RstAquifer(std::shared_ptr rstView, + const EclipseGrid* grid, + const UnitSystem& usys); + + RstAquifer(const RstAquifer& rhs); + RstAquifer(RstAquifer&& rhs); + RstAquifer& operator=(const RstAquifer& rhs); + RstAquifer& operator=(RstAquifer&& rhs); + + ~RstAquifer(); + + bool hasAnalyticAquifers() const; + + const std::vector& carterTracy() const; + const std::vector& fetkovich() const; + const std::unordered_map& connections() const; + + private: + class Implementation; + std::unique_ptr pImpl_; + }; +}} // Opm::RestartIO + +#endif // OPM_RESTART_AQUIFER_HPP diff --git a/opm/io/eclipse/rst/state.hpp b/opm/io/eclipse/rst/state.hpp index a9767f074..288e64ad5 100644 --- a/opm/io/eclipse/rst/state.hpp +++ b/opm/io/eclipse/rst/state.hpp @@ -24,6 +24,7 @@ #include #include +#include #include #include #include @@ -32,6 +33,9 @@ #include +namespace Opm { + class EclipseGrid; +} // namespace Opm namespace Opm { namespace EclIO { class RestartFileView; @@ -40,17 +44,17 @@ namespace Opm { namespace EclIO { namespace Opm { namespace RestartIO { struct RstState { - RstState(const ::Opm::UnitSystem& unit_system, - const std::vector& intehead, - const std::vector& logihead, - const std::vector& doubhead); + RstState(std::shared_ptr rstView, + const ::Opm::EclipseGrid* grid); - static RstState load(std::shared_ptr rstView); + static RstState load(std::shared_ptr rstView, + const ::Opm::EclipseGrid* grid = nullptr); const RstWell& get_well(const std::string& wname) const; - const ::Opm::UnitSystem unit_system; + ::Opm::UnitSystem unit_system; RstHeader header; + RstAquifer aquifers; std::vector wells; std::vector groups; std::vector udqs; diff --git a/opm/parser/eclipse/EclipseState/Aquifer/Aquancon.hpp b/opm/parser/eclipse/EclipseState/Aquifer/Aquancon.hpp index 628cf66e6..8eb75203d 100644 --- a/opm/parser/eclipse/EclipseState/Aquifer/Aquancon.hpp +++ b/opm/parser/eclipse/EclipseState/Aquifer/Aquancon.hpp @@ -37,6 +37,10 @@ namespace Opm { class Deck; } +namespace Opm { namespace RestartIO { + class RstAquifer; +}} // Opm::RestartIO + namespace Opm { class Aquancon { @@ -82,11 +86,12 @@ namespace Opm { } }; - Aquancon() = default; Aquancon(const EclipseGrid& grid, const Deck& deck); explicit Aquancon(const std::unordered_map>& data); + void loadFromRestart(const RestartIO::RstAquifer& rst_aquifers); + static Aquancon serializeObject(); const std::unordered_map>& data() const; diff --git a/opm/parser/eclipse/EclipseState/Aquifer/AquiferCT.hpp b/opm/parser/eclipse/EclipseState/Aquifer/AquiferCT.hpp index 2f8f89072..82792734a 100644 --- a/opm/parser/eclipse/EclipseState/Aquifer/AquiferCT.hpp +++ b/opm/parser/eclipse/EclipseState/Aquifer/AquiferCT.hpp @@ -36,6 +36,10 @@ namespace Opm { class TableManager; } +namespace Opm { namespace RestartIO { + class RstAquifer; +}} // Opm::RestartIO + namespace Opm { class AquiferCT { @@ -43,6 +47,8 @@ namespace Opm { struct AQUCT_data { + friend class AquiferCT; + AQUCT_data() = default; AQUCT_data(const DeckRecord& record, const TableManager& tables); AQUCT_data(const int aqID, @@ -117,6 +123,9 @@ namespace Opm { AquiferCT(const TableManager& tables, const Deck& deck); AquiferCT(const std::vector& data); + void loadFromRestart(const RestartIO::RstAquifer& rst, + const TableManager& tables); + static AquiferCT serializeObject(); std::size_t size() const; diff --git a/opm/parser/eclipse/EclipseState/Aquifer/AquiferConfig.hpp b/opm/parser/eclipse/EclipseState/Aquifer/AquiferConfig.hpp index 0ee07de33..bfdc9feb1 100644 --- a/opm/parser/eclipse/EclipseState/Aquifer/AquiferConfig.hpp +++ b/opm/parser/eclipse/EclipseState/Aquifer/AquiferConfig.hpp @@ -26,11 +26,17 @@ #include namespace Opm { + class TableManager; + class EclipseGrid; + class Deck; + class FieldPropsManager; +} // namespace Opm -class TableManager; -class EclipseGrid; -class Deck; -class FieldPropsManager; +namespace Opm { namespace RestartIO { + class RstAquifer; +}} // namespace Opm::RestartIO + +namespace Opm { class AquiferConfig { public: @@ -40,6 +46,8 @@ public: AquiferConfig(const Aquifetp& fetp, const AquiferCT& ct, const Aquancon& conn); void load_connections(const Deck& deck, const EclipseGrid& grid); + void loadFromRestart(const RestartIO::RstAquifer& aquifers, + const TableManager& tables); static AquiferConfig serializeObject(); diff --git a/opm/parser/eclipse/EclipseState/Aquifer/Aquifetp.hpp b/opm/parser/eclipse/EclipseState/Aquifer/Aquifetp.hpp index 5e165e8f5..b280a6e58 100644 --- a/opm/parser/eclipse/EclipseState/Aquifer/Aquifetp.hpp +++ b/opm/parser/eclipse/EclipseState/Aquifer/Aquifetp.hpp @@ -36,6 +36,10 @@ namespace Opm { class TableManager; } +namespace Opm { namespace RestartIO { + class RstAquifer; +}} // Opm::RestartIO + namespace Opm { class Aquifetp { @@ -43,6 +47,8 @@ class Aquifetp { struct AQUFETP_data { + friend class Aquifetp; + AQUFETP_data() = default; AQUFETP_data(const DeckRecord& record, const TableManager& tables); AQUFETP_data(const int aquiferID_, @@ -98,6 +104,9 @@ class Aquifetp { Aquifetp(const TableManager& tables, const Deck& deck); explicit Aquifetp(const std::vector& data); + void loadFromRestart(const RestartIO::RstAquifer& rst, + const TableManager& tables); + static Aquifetp serializeObject(); const std::vector& data() const; diff --git a/opm/parser/eclipse/EclipseState/EclipseState.hpp b/opm/parser/eclipse/EclipseState/EclipseState.hpp index 16038f26b..32c331162 100644 --- a/opm/parser/eclipse/EclipseState/EclipseState.hpp +++ b/opm/parser/eclipse/EclipseState/EclipseState.hpp @@ -15,7 +15,7 @@ You should have received a copy of the GNU General Public License along with OPM. If not, see . - */ +*/ #ifndef OPM_ECLIPSE_STATE_HPP #define OPM_ECLIPSE_STATE_HPP @@ -37,13 +37,18 @@ #include namespace Opm { - class Deck; class DeckKeyword; class InitConfig; class IOConfig; class DeckSection; +} // namespace Opm +namespace Opm { namespace RestartIO { + class RstAquifer; +}} // namespace Opm::RestartIO + +namespace Opm { class EclipseState { public: @@ -100,6 +105,8 @@ namespace Opm { const AquiferConfig& aquifer() const; const TracerConfig& tracer() const; + void loadRestartAquifers(const RestartIO::RstAquifer& aquifers); + template void serializeOp(Serializer& serializer) { @@ -148,6 +155,6 @@ namespace Opm { std::string m_title; TracerConfig tracer_config; }; -} +} // namespace Opm #endif // OPM_ECLIPSE_STATE_HPP diff --git a/src/opm/io/eclipse/RestartFileView.cpp b/src/opm/io/eclipse/RestartFileView.cpp index 985d92a4c..c0b05892c 100644 --- a/src/opm/io/eclipse/RestartFileView.cpp +++ b/src/opm/io/eclipse/RestartFileView.cpp @@ -132,6 +132,32 @@ public: return this->getKeyword(ihkw, 0); } + const std::vector& logihead() + { + const auto lhkw = std::string { "LOGIHEAD" }; + + if (! this->hasKeyword(lhkw)) { + throw std::domain_error { + "Purported Restart File Does not Have Logical Header" + }; + } + + return this->getKeyword(lhkw, 0); + } + + const std::vector& doubhead() + { + const auto dhkw = std::string { "DOUBHEAD" }; + + if (! this->hasKeyword(dhkw)) { + throw std::domain_error { + "Purported Restart File Does not Have Double Header" + }; + } + + return this->getKeyword(dhkw, 0); + } + private: using RstFile = std::shared_ptr; @@ -170,7 +196,6 @@ Implementation(std::shared_ptr restart_file, const auto& type = std::get<1>(vector); switch (type) { - case ::Opm::EclIO::eclArrType::LOGI: case ::Opm::EclIO::eclArrType::MESS: // Currently ignored continue; @@ -240,6 +265,16 @@ const std::vector& Opm::EclIO::RestartFileView::intehead() const return this->pImpl_->intehead(); } +const std::vector& Opm::EclIO::RestartFileView::logihead() const +{ + return this->pImpl_->logihead(); +} + +const std::vector& Opm::EclIO::RestartFileView::doubhead() const +{ + return this->pImpl_->doubhead(); +} + template bool Opm::EclIO::RestartFileView::hasKeyword(const std::string& vector) const { diff --git a/src/opm/io/eclipse/rst/aquifer.cpp b/src/opm/io/eclipse/rst/aquifer.cpp new file mode 100644 index 000000000..0d4fc71f6 --- /dev/null +++ b/src/opm/io/eclipse/rst/aquifer.cpp @@ -0,0 +1,683 @@ +/* + Copyright 2021 Equinor ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#include + +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include + +namespace VI = ::Opm::RestartIO::Helpers::VectorItems; + +namespace { +template +boost::iterator_range::const_iterator> +getDataWindow(const std::vector& arr, + const std::size_t windowSize, + const std::size_t entity, + const std::size_t subEntity = 0, + const std::size_t maxSubEntitiesPerEntity = 1) +{ + const auto off = + windowSize * (subEntity + maxSubEntitiesPerEntity*entity); + + auto begin = arr.begin() + off; + auto end = begin + windowSize; + + return { begin, end }; +} + +// --------------------------------------------------------------------- + +class ConnectionVectors +{ +public: + template + using Window = boost::iterator_range< + typename std::vector::const_iterator + >; + + explicit ConnectionVectors(const std::vector& intehead, + std::shared_ptr rst_view); + + Window icaq(const int occurrence, + const std::size_t connectionID) const; + + Window scaq(const int occurrence, + const std::size_t connectionID) const; + +private: + std::size_t numIConnElem_; + std::size_t numSConnElem_; + + std::shared_ptr rstView_; +}; + +ConnectionVectors::ConnectionVectors(const std::vector& intehead, + std::shared_ptr rst_view) + : numIConnElem_(intehead[VI::intehead::NICAQZ]) + , numSConnElem_(intehead[VI::intehead::NSCAQZ]) + , rstView_ (std::move(rst_view)) +{} + +ConnectionVectors::Window +ConnectionVectors::icaq(const int occurrence, + const std::size_t connectionID) const +{ + return getDataWindow(this->rstView_->getKeyword("ICAQ", occurrence), + this->numIConnElem_, connectionID); +} + +ConnectionVectors::Window +ConnectionVectors::scaq(const int occurrence, + const std::size_t connectionID) const +{ + return getDataWindow(this->rstView_->getKeyword("SCAQ", occurrence), + this->numSConnElem_, connectionID); +} + +// --------------------------------------------------------------------- + +class AquiferVectors +{ +public: + template + using Window = boost::iterator_range< + typename std::vector::const_iterator + >; + + explicit AquiferVectors(const std::vector& intehead, + std::shared_ptr rst_view); + + int maxAquiferID() const; + + Window iaaq(const std::size_t aquiferID) const; + Window saaq(const std::size_t aquiferID) const; + Window xaaq(const std::size_t aquiferID) const; + +private: + int maxAquiferID_{}; + std::size_t numIAquifElem_; + std::size_t numSAquifElem_; + std::size_t numXAquifElem_; + + std::shared_ptr rstView_; +}; + +AquiferVectors::AquiferVectors(const std::vector& intehead, + std::shared_ptr rst_view) + : maxAquiferID_ (intehead[VI::intehead::MAX_AN_AQUIFER_ID]) + , numIAquifElem_(intehead[VI::intehead::NIAAQZ]) + , numSAquifElem_(intehead[VI::intehead::NSAAQZ]) + , numXAquifElem_(intehead[VI::intehead::NXAAQZ]) + , rstView_ (std::move(rst_view)) +{} + +int AquiferVectors::maxAquiferID() const +{ + return this->maxAquiferID_; +} + +AquiferVectors::Window +AquiferVectors::iaaq(const std::size_t aquiferID) const +{ + return getDataWindow(this->rstView_->getKeyword("IAAQ"), + this->numIAquifElem_, aquiferID); +} + +AquiferVectors::Window +AquiferVectors::saaq(const std::size_t aquiferID) const +{ + return getDataWindow(this->rstView_->getKeyword("SAAQ"), + this->numSAquifElem_, aquiferID); +} + +AquiferVectors::Window +AquiferVectors::xaaq(const std::size_t aquiferID) const +{ + return getDataWindow(this->rstView_->getKeyword("XAAQ"), + this->numXAquifElem_, aquiferID); +} + +// --------------------------------------------------------------------- + +class ConnectionOccurrence +{ +public: + explicit ConnectionOccurrence(const int numAquifers, + Opm::EclIO::RestartFileView& rst_view); + + int icaq(const int aquiferID) const + { + return this->occurrence_[aquiferID].I; + } + + int scaq(const int aquiferID) const + { + return this->occurrence_[aquiferID].S; + } + +private: + struct Occurrence + { + int I{-1}; + int S{-1}; + }; + + std::vector occurrence_{}; +}; + +ConnectionOccurrence::ConnectionOccurrence(const int maxAquiferID, + Opm::EclIO::RestartFileView& rst_view) + : occurrence_(maxAquiferID) +{ + const auto integerAquifID = std::string{ "ICAQNUM" }; + const auto floatAquifID = std::string{ "SCAQNUM" }; + + if ((rst_view.occurrenceCount(integerAquifID) < maxAquiferID) || + (rst_view.occurrenceCount(floatAquifID) < maxAquiferID)) + { + throw std::invalid_argument { + "Structural inconsistency for analytic " + "aquifer connections in restart file" + }; + } + + // Note: + // One separate xCAQNUM table for each aquifer ID in 1..maxAquiferID. + // + // On the one hand, Flow's output system at the time of writing + // guarantees that + // + // occurrence_[i].I == occurrence_[i].S == i + // + // for all 'i' in 0..maxAquiferID-1. On the other hand, the file + // format allows for general ordering. This mapping facility exists + // mainly to handle the general case. + for (auto occurrence = 0*maxAquiferID; occurrence < maxAquiferID; ++occurrence) { + const auto& iAquifID = rst_view.getKeyword(integerAquifID, occurrence); + const auto& sAquifID = rst_view.getKeyword(floatAquifID, occurrence); + + this->occurrence_[iAquifID.front() - 1].I = occurrence; + this->occurrence_[sAquifID.front() - 1].S = occurrence; + } +} + +// --------------------------------------------------------------------- + +Opm::FaceDir::DirEnum face_direction(const int directionValue) +{ + using FDValue = VI::IAnalyticAquiferConn::Value::FaceDirection; + using FDir = Opm::FaceDir::DirEnum; + + switch (directionValue) { + case FDValue::IMinus: return FDir::XMinus; + case FDValue::IPlus: return FDir::XPlus; + case FDValue::JMinus: return FDir::YMinus; + case FDValue::JPlus: return FDir::YPlus; + case FDValue::KMinus: return FDir::ZMinus; + case FDValue::KPlus: return FDir::ZPlus; + } + + throw std::invalid_argument { + fmt::format("Unknown Face Direction {}", directionValue) + }; +} + +template +Opm::FaceDir::DirEnum face_direction(const IcaqArray& icaq) +{ + using Ix = VI::IAnalyticAquiferConn::index; + + return face_direction(icaq[Ix::FaceDirection]); +} + +template +std::size_t identify_global_cell(const Opm::EclipseGrid& grid, + const IcaqArray& icaq) +{ + using Ix = VI::IAnalyticAquiferConn::index; + + const auto i = static_cast(icaq[Ix::Index_I] - 1); + const auto j = static_cast(icaq[Ix::Index_J] - 1); + const auto k = static_cast(icaq[Ix::Index_K] - 1); + + return grid.getGlobalIndex(i, j, k); +} + +template +double influx_coefficient(const ScaqArray scaq, const double tot_influx) +{ + using Ix = VI::SAnalyticAquiferConn::index; + + return scaq[Ix::InfluxFraction] * tot_influx; +} + +template +double effective_facearea(const ScaqArray scaq, const double tot_influx) +{ + using Ix = VI::SAnalyticAquiferConn::index; + + return scaq[Ix::FaceAreaToInfluxCoeff] * tot_influx; +} + +template +void load_analytic_aquifer_cell(const IcaqArray icaq, + const ScaqArray scaq, + const Opm::EclipseGrid& grid, + const double tot_influx, + Opm::RestartIO::RstAquifer::Connections& connections) +{ + const auto glob_cell = identify_global_cell(grid, icaq); + const auto influx_coeff = influx_coefficient(scaq, tot_influx); + const auto eff_facearea = effective_facearea(scaq, tot_influx); + const auto face_dir = face_direction(icaq); + + connections.emplace_back(glob_cell, influx_coeff, eff_facearea, face_dir); +} + +Opm::RestartIO::RstAquifer::Connections +load_analytic_aquifer_cells(const ConnectionOccurrence& occurence, + const ConnectionVectors& connections, + const Opm::EclipseGrid& grid, + const std::size_t num_conn, + const double tot_influx, + const int aquifer_id) +{ + auto cells = Opm::RestartIO::RstAquifer::Connections{}; + cells.reserve(num_conn); + + const auto occur_icaq = occurence.icaq(aquifer_id); + const auto occur_scaq = occurence.scaq(aquifer_id); + + for (auto connID = 0*num_conn; connID < num_conn; ++connID) { + auto icaq = connections.icaq(occur_icaq, connID); + auto scaq = connections.scaq(occur_scaq, connID); + + load_analytic_aquifer_cell(icaq, scaq, grid, tot_influx, cells); + } + + return cells; +} + +// --------------------------------------------------------------------- + +Opm::RestartIO::RstAquifer::CarterTracy +load_carter_tracy(const int aquiferID, + const AquiferVectors& aquifers, + const Opm::UnitSystem& usys) +{ + auto aquifer = Opm::RestartIO::RstAquifer::CarterTracy{}; + + using M = Opm::UnitSystem::measure; + using IntIX = VI::IAnalyticAquifer::index; + using RealIX = VI::SAnalyticAquifer::index; + using DoubIX = VI::XAnalyticAquifer::index; + + const auto iaaq = aquifers.iaaq(aquiferID); + const auto saaq = aquifers.saaq(aquiferID); + const auto xaaq = aquifers.xaaq(aquiferID); + + aquifer.aquiferID = aquiferID + 1; + aquifer.inftableID = iaaq[IntIX::CTInfluenceFunction]; + aquifer.pvttableID = iaaq[IntIX::WatPropTable]; + + aquifer.porosity = usys.to_si(M::identity, saaq[RealIX::CTPorosity]); + aquifer.datum_depth = usys.to_si(M::length, saaq[RealIX::DatumDepth]); + + // Note: *from_si()* to work around the fact that we don't have a + // compressibility unit. + aquifer.total_compr = + usys.from_si(M::pressure, saaq[RealIX::Compressibility]); + + aquifer.permeability = + usys.to_si(M::permeability, saaq[RealIX::CTPermeability]); + + aquifer.inner_radius = usys.to_si(M::length, saaq[RealIX::CTRadius]); + aquifer.thickness = usys.to_si(M::length, saaq[RealIX::CTThickness]); + aquifer.angle_fraction = usys.to_si(M::identity, saaq[RealIX::CTAngle]); + aquifer.initial_pressure = usys.to_si(M::pressure, saaq[RealIX::InitPressure]); + aquifer.time_constant = usys.to_si(M::time, 1.0 / xaaq[DoubIX::CTRecipTimeConst]); + + // Note: *from_si()* for the pressure unit here since 'beta' is total + // influx (volume) per unit pressure drop. + aquifer.influx_constant = + usys.to_si(M::volume, usys.from_si(M::pressure, xaaq[DoubIX::CTInfluxConstant])); + + aquifer.water_density = usys.to_si(M::density, saaq[RealIX::CTWatMassDensity]); + aquifer.water_viscosity = usys.to_si(M::viscosity, saaq[RealIX::CTWatViscosity]); + + return aquifer; +} + +// --------------------------------------------------------------------- + +Opm::RestartIO::RstAquifer::Fetkovich +load_fetkovich(const int aquiferID, + const AquiferVectors& aquifers, + const Opm::UnitSystem& usys) +{ + auto aquifer = Opm::RestartIO::RstAquifer::Fetkovich{}; + + using M = Opm::UnitSystem::measure; + using IntIX = VI::IAnalyticAquifer::index; + using RealIX = VI::SAnalyticAquifer::index; + + const auto iaaq = aquifers.iaaq(aquiferID); + const auto saaq = aquifers.saaq(aquiferID); + + aquifer.aquiferID = aquiferID + 1; + aquifer.pvttableID = iaaq[IntIX::WatPropTable]; + + aquifer.prod_index = + usys.to_si(M::liquid_productivity_index, saaq[RealIX::FetProdIndex]); + + // Note: *from_si()* to work around the fact that we don't have a + // compressibility unit. + aquifer.total_compr = + usys.from_si(M::pressure, saaq[RealIX::Compressibility]); + + aquifer.initial_watvolume = + usys.to_si(M::liquid_surface_volume, saaq[RealIX::FetInitVol]); + + aquifer.datum_depth = usys.to_si(M::length, saaq[RealIX::DatumDepth]); + aquifer.initial_pressure = usys.to_si(M::pressure, saaq[RealIX::InitPressure]); + aquifer.time_constant = usys.to_si(M::time, saaq[RealIX::FetTimeConstant]); + + return aquifer; +} + +// --------------------------------------------------------------------- + +int num_aquifers(const std::vector& intehead) +{ + return intehead[VI::intehead::NAQUIF]; +} + +int num_aquifer_connections(const AquiferVectors& aquifers, + const std::size_t aquiferID) +{ + using Ix = VI::IAnalyticAquifer::index; + + auto iaaq = aquifers.iaaq(aquiferID); + + return iaaq[Ix::NumAquiferConn]; +} + +std::unordered_map +load_aquifer_connections(const ConnectionOccurrence& occurence, + const AquiferVectors& aquifers, + const ConnectionVectors& connections, + const Opm::EclipseGrid& grid, + const Opm::UnitSystem& usys, + const int max_aquifer_id) +{ + auto aqConn = std::unordered_map{}; + + auto tot_influx = [&aquifers, &usys](const std::size_t aquiferID) -> double + { + using M = Opm::UnitSystem::measure; + using Ix = VI::XAnalyticAquifer::index; + + auto xaaq = aquifers.xaaq(aquiferID); + + return usys.to_si(M::length, usys.to_si(M::length, xaaq[Ix::TotalInfluxCoeff])); + }; + + auto load = [&occurence, &connections, &grid] + (const std::size_t num_connections, + const double total_influx, + const int aquifer_id) + { + return load_analytic_aquifer_cells(occurence, connections, grid, + num_connections, total_influx, aquifer_id); + }; + + for (auto aquiferID = 0*max_aquifer_id; aquiferID < max_aquifer_id; ++aquiferID) { + const auto num_connections = num_aquifer_connections(aquifers, aquiferID); + + if (num_connections == 0) { + continue; + } + + const auto total_influx = tot_influx(aquiferID); + aqConn.insert_or_assign(aquiferID + 1, load(num_connections, total_influx, aquiferID)); + } + + return aqConn; +} +} // Anonymous + +class Opm::RestartIO::RstAquifer::Implementation +{ +public: + explicit Implementation(std::shared_ptr rstView, + const EclipseGrid* grid, + const UnitSystem& usys); + + bool hasAnalyticAquifers() const + { + return ! (this->connections_.empty() && + this->carterTracy_.empty() && + this->fetkovich_ .empty()); + } + + const std::vector& carterTracy() const + { + return this->carterTracy_; + } + + const std::vector& fetkovich() const + { + return this->fetkovich_; + } + + const std::unordered_map& connections() const + { + return this->connections_; + } + +private: + std::unordered_map connections_{}; + std::vector carterTracy_{}; + std::vector fetkovich_{}; + + void loadAnalyticAquiferConnections(const AquiferVectors& vectors, + std::shared_ptr rstView, + const EclipseGrid& grid, + const UnitSystem& usys); + + void loadAnalyticAquifers(const AquiferVectors& aquifers, + const UnitSystem& usys); + void loadAnalyticAquifer(const int aquiferID, + const AquiferVectors& aquifers, + const UnitSystem& usys); + + void loadCarterTracy(const int aquiferID, + const AquiferVectors& aquifers, + const UnitSystem& usys); + + void loadFetkovich(const int aquiferID, + const AquiferVectors& aquifers, + const UnitSystem& usys); +}; + +Opm::RestartIO::RstAquifer::Implementation:: +Implementation(std::shared_ptr rstView, + const EclipseGrid* grid, + const UnitSystem& usys) +{ + const auto numAquifers = num_aquifers(rstView->intehead()); + if ((numAquifers == 0) || (grid == nullptr)) { + return; + } + + const auto aquifers = AquiferVectors { rstView->intehead(), rstView }; + + this->loadAnalyticAquiferConnections(aquifers, std::move(rstView), *grid, usys); + this->loadAnalyticAquifers(aquifers, usys); +} + +void +Opm::RestartIO::RstAquifer::Implementation:: +loadAnalyticAquiferConnections(const AquiferVectors& aquifers, + std::shared_ptr rstView, + const EclipseGrid& grid, + const UnitSystem& usys) +{ + const auto occurence = ConnectionOccurrence { aquifers.maxAquiferID(), *rstView }; + const auto connections = ConnectionVectors { rstView->intehead(), std::move(rstView) }; + + this->connections_ = + load_aquifer_connections(occurence, aquifers, connections, + grid, usys, aquifers.maxAquiferID()); +} + +void +Opm::RestartIO::RstAquifer::Implementation:: +loadAnalyticAquifers(const AquiferVectors& aquifers, + const UnitSystem& usys) +{ + const auto maxAquiferID = aquifers.maxAquiferID(); + for (auto aquiferID = 0*maxAquiferID; aquiferID < maxAquiferID; ++aquiferID) { + if (num_aquifer_connections(aquifers, aquiferID) == 0) { + // Skip aquifers without connections. Likely to be sparsely + // allocated aquifer IDs (e.g., 1, 2, 5). + continue; + } + + this->loadAnalyticAquifer(aquiferID, aquifers, usys); + } +} + +void +Opm::RestartIO::RstAquifer::Implementation:: +loadAnalyticAquifer(const int aquiferID, + const AquiferVectors& aquifers, + const UnitSystem& usys) +{ + const auto iaaq = aquifers.iaaq(aquiferID); + const auto type = iaaq[VI::IAnalyticAquifer::index::TypeRelated1]; + + switch (type) { + case VI::IAnalyticAquifer::Value::ModelType::CarterTracy: + this->loadCarterTracy(aquiferID, aquifers, usys); + return; + + case VI::IAnalyticAquifer::Value::ModelType::Fetkovich: + this->loadFetkovich(aquiferID, aquifers, usys); + return; + } + + throw std::invalid_argument { + fmt::format("Analytic aquifer {} (type {}) is neither Carter-Tracy " + "nor Fetkovich in restart input file", aquiferID + 1, type) + }; +} + +void +Opm::RestartIO::RstAquifer::Implementation:: +loadCarterTracy(const int aquiferID, + const AquiferVectors& aquifers, + const UnitSystem& usys) +{ + this->carterTracy_.push_back(load_carter_tracy(aquiferID, aquifers, usys)); +} + +void +Opm::RestartIO::RstAquifer::Implementation:: +loadFetkovich(const int aquiferID, + const AquiferVectors& aquifers, + const UnitSystem& usys) +{ + this->fetkovich_.push_back(load_fetkovich(aquiferID, aquifers, usys)); +} + +// --------------------------------------------------------------------- + +Opm::RestartIO::RstAquifer::RstAquifer(std::shared_ptr rstView, + const EclipseGrid* grid, + const UnitSystem& usys) + : pImpl_{ new Implementation{ std::move(rstView), grid, usys } } +{} + +Opm::RestartIO::RstAquifer::RstAquifer(const RstAquifer& rhs) + : pImpl_{ new Implementation{ *rhs.pImpl_ } } +{} + +Opm::RestartIO::RstAquifer::RstAquifer(RstAquifer&& rhs) + : pImpl_{ std::move(rhs.pImpl_) } +{} + +Opm::RestartIO::RstAquifer& +Opm::RestartIO::RstAquifer::operator=(const RstAquifer& rhs) +{ + this->pImpl_.reset(new Implementation{ *rhs.pImpl_ }); + return *this; +} + +Opm::RestartIO::RstAquifer& +Opm::RestartIO::RstAquifer::operator=(RstAquifer&& rhs) +{ + this->pImpl_ = std::move(rhs.pImpl_); + return *this; +} + +Opm::RestartIO::RstAquifer::~RstAquifer() +{} + +bool Opm::RestartIO::RstAquifer::hasAnalyticAquifers() const +{ + return this->pImpl_->hasAnalyticAquifers(); +} + +const std::vector& +Opm::RestartIO::RstAquifer::carterTracy() const +{ + return this->pImpl_->carterTracy(); +} + +const std::vector& +Opm::RestartIO::RstAquifer::fetkovich() const +{ + return this->pImpl_->fetkovich(); +} + +const std::unordered_map& +Opm::RestartIO::RstAquifer::connections() const +{ + return this->pImpl_->connections(); +} diff --git a/src/opm/io/eclipse/rst/state.cpp b/src/opm/io/eclipse/rst/state.cpp index 37a958356..d7153fa7c 100644 --- a/src/opm/io/eclipse/rst/state.cpp +++ b/src/opm/io/eclipse/rst/state.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include @@ -64,14 +65,13 @@ namespace VI = ::Opm::RestartIO::Helpers::VectorItems; namespace Opm { namespace RestartIO { -RstState::RstState(const ::Opm::UnitSystem& unit_system_, - const std::vector& intehead, - const std::vector& logihead, - const std::vector& doubhead): - unit_system(unit_system_), - header(unit_system_, intehead, logihead, doubhead) +RstState::RstState(std::shared_ptr rstView, + const ::Opm::EclipseGrid* grid) + : unit_system(rstView->intehead()[VI::intehead::UNIT]) + , header(unit_system, rstView->intehead(), rstView->logihead(), rstView->doubhead()) + , aquifers(rstView, grid, unit_system) { - this->load_tuning(intehead, doubhead); + this->load_tuning(rstView->intehead(), rstView->doubhead()); } @@ -273,15 +273,11 @@ const RstWell& RstState::get_well(const std::string& wname) const { return *well_iter; } -RstState RstState::load(std::shared_ptr rstView) { - const auto& intehead = rstView->getKeyword("INTEHEAD"); - const auto& logihead = rstView->getKeyword("LOGIHEAD"); - const auto& doubhead = rstView->getKeyword("DOUBHEAD"); +RstState RstState::load(std::shared_ptr rstView, + const ::Opm::EclipseGrid* grid) +{ + RstState state(rstView, grid); - auto unit_id = intehead[VI::intehead::UNIT]; - ::Opm::UnitSystem unit_system(unit_id); - - RstState state(unit_system, intehead, logihead, doubhead); if (state.header.ngroup > 0) { const auto& zgrp = rstView->getKeyword("ZGRP"); const auto& igrp = rstView->getKeyword("IGRP"); diff --git a/src/opm/parser/eclipse/EclipseState/Aquifer/Aquancon.cpp b/src/opm/parser/eclipse/EclipseState/Aquifer/Aquancon.cpp index b4ae3ea57..5a484ed24 100644 --- a/src/opm/parser/eclipse/EclipseState/Aquifer/Aquancon.cpp +++ b/src/opm/parser/eclipse/EclipseState/Aquifer/Aquancon.cpp @@ -15,10 +15,12 @@ You should have received a copy of the GNU General Public License along with OPM. If not, see . - */ +*/ #include +#include + #include #include @@ -46,6 +48,36 @@ #include "AquiferHelpers.hpp" +namespace { + Opm::Aquancon::AquancCell + makeAquiferCell(const int aquiferID, + const Opm::RestartIO::RstAquifer::Connections::Cell& rst_cell) + { + return { + aquiferID, + rst_cell.global_index, + rst_cell.influx_coeff, + rst_cell.effective_facearea, + rst_cell.face_dir + }; + } + + std::vector + makeAquiferConnections(const int aquiferID, + const Opm::RestartIO::RstAquifer::Connections& rst_connections) + { + auto connections = std::vector{}; + const auto& rst_cells = rst_connections.cells(); + + connections.reserve(rst_cells.size()); + for (const auto& rst_cell : rst_cells) { + connections.push_back(makeAquiferCell(aquiferID, rst_cell)); + } + + return connections; + } +} // Anonymous + namespace Opm { namespace{ @@ -183,6 +215,15 @@ namespace Opm { cells(data) {} + void Aquancon::loadFromRestart(const RestartIO::RstAquifer& rst_aquifers) + { + this->cells.clear(); + + for (const auto& [aquiferID, rst_connections] : rst_aquifers.connections()) { + this->cells.insert_or_assign(aquiferID, makeAquiferConnections(aquiferID, rst_connections)); + } + } + const std::unordered_map>& Aquancon::data() const { return this->cells; } diff --git a/src/opm/parser/eclipse/EclipseState/Aquifer/AquiferCT.cpp b/src/opm/parser/eclipse/EclipseState/Aquifer/AquiferCT.cpp index 43172d506..b3c763405 100644 --- a/src/opm/parser/eclipse/EclipseState/Aquifer/AquiferCT.cpp +++ b/src/opm/parser/eclipse/EclipseState/Aquifer/AquiferCT.cpp @@ -15,10 +15,12 @@ You should have received a copy of the GNU General Public License along with OPM. If not, see . - */ +*/ #include +#include + #include #include @@ -42,6 +44,29 @@ #include #include +namespace { + Opm::AquiferCT::AQUCT_data + makeAquifer(const Opm::RestartIO::RstAquifer::CarterTracy& rst_aquifer) + { + auto aquifer = Opm::AquiferCT::AQUCT_data{}; + + aquifer.aquiferID = rst_aquifer.aquiferID; + aquifer.inftableID = rst_aquifer.inftableID; + aquifer.pvttableID = rst_aquifer.pvttableID; + + aquifer.porosity = rst_aquifer.porosity; + aquifer.datum_depth = rst_aquifer.datum_depth; + aquifer.total_compr = rst_aquifer.total_compr; + aquifer.inner_radius = rst_aquifer.inner_radius; + aquifer.permeability = rst_aquifer.permeability; + aquifer.thickness = rst_aquifer.thickness; + aquifer.angle_fraction = rst_aquifer.angle_fraction; + aquifer.initial_pressure = rst_aquifer.initial_pressure; + + return aquifer; + } +} + namespace Opm { namespace { @@ -199,6 +224,30 @@ AquiferCT::AquiferCT(const std::vector& data) : m_aquct(data) {} +void AquiferCT::loadFromRestart(const RestartIO::RstAquifer& rst, + const TableManager& tables) +{ + this->m_aquct.clear(); + + const auto& rst_aquifers = rst.carterTracy(); + if (! rst_aquifers.empty()) { + this->m_aquct.reserve(rst_aquifers.size()); + } + + for (const auto& rst_aquifer : rst_aquifers) { + this->m_aquct.push_back(makeAquifer(rst_aquifer)); + + auto& new_aquifer = this->m_aquct.back(); + new_aquifer.finishInitialisation(tables); + + // Assign 'private:' members in AQUCT_data. + new_aquifer.time_constant_ = rst_aquifer.time_constant; + new_aquifer.influx_constant_ = rst_aquifer.influx_constant; + new_aquifer.water_density_ = rst_aquifer.water_density; + new_aquifer.water_viscosity_ = rst_aquifer.water_viscosity; + } +} + AquiferCT AquiferCT::serializeObject() { return { { AQUCT_data::serializeObject() } }; diff --git a/src/opm/parser/eclipse/EclipseState/Aquifer/AquiferConfig.cpp b/src/opm/parser/eclipse/EclipseState/Aquifer/AquiferConfig.cpp index 02c567724..a7b61396a 100644 --- a/src/opm/parser/eclipse/EclipseState/Aquifer/AquiferConfig.cpp +++ b/src/opm/parser/eclipse/EclipseState/Aquifer/AquiferConfig.cpp @@ -47,6 +47,14 @@ void AquiferConfig::load_connections(const Deck& deck, const EclipseGrid& grid) this->aqconn = Aquancon(grid, deck); } +void AquiferConfig::loadFromRestart(const RestartIO::RstAquifer& aquifers, + const TableManager& tables) +{ + this->aquifetp.loadFromRestart(aquifers, tables); + this->aquiferct.loadFromRestart(aquifers, tables); + this->aqconn.loadFromRestart(aquifers); +} + AquiferConfig AquiferConfig::serializeObject() { AquiferConfig result; diff --git a/src/opm/parser/eclipse/EclipseState/Aquifer/Aquifetp.cpp b/src/opm/parser/eclipse/EclipseState/Aquifer/Aquifetp.cpp index bd96c798c..c90079678 100644 --- a/src/opm/parser/eclipse/EclipseState/Aquifer/Aquifetp.cpp +++ b/src/opm/parser/eclipse/EclipseState/Aquifer/Aquifetp.cpp @@ -19,6 +19,8 @@ #include +#include + #include #include @@ -32,6 +34,25 @@ #include +namespace { + Opm::Aquifetp::AQUFETP_data + makeAquifer(const Opm::RestartIO::RstAquifer::Fetkovich& rst_aquifer) + { + auto aquifer = Opm::Aquifetp::AQUFETP_data{}; + + aquifer.aquiferID = rst_aquifer.aquiferID; + aquifer.pvttableID = rst_aquifer.pvttableID; + + aquifer.prod_index = rst_aquifer.prod_index; + aquifer.total_compr = rst_aquifer.total_compr; + aquifer.initial_watvolume = rst_aquifer.initial_watvolume; + aquifer.datum_depth = rst_aquifer.datum_depth; + aquifer.initial_pressure = rst_aquifer.initial_pressure; + + return aquifer; + } +} + namespace Opm { using AQUFETP = ParserKeywords::AQUFETP; @@ -131,6 +152,26 @@ Aquifetp::Aquifetp(const std::vector& data) : m_aqufetp(data) {} +void Aquifetp::loadFromRestart(const RestartIO::RstAquifer& rst, + const TableManager& tables) +{ + this->m_aqufetp.clear(); + + const auto& rst_aquifers = rst.fetkovich(); + if (! rst_aquifers.empty()) { + this->m_aqufetp.reserve(rst_aquifers.size()); + } + + for (const auto& rst_aquifer : rst_aquifers) { + this->m_aqufetp.push_back(makeAquifer(rst_aquifer)); + + auto& new_aquifer = this->m_aqufetp.back(); + new_aquifer.finishInitialisation(tables); + + // Assign 'private:' members in AQUFETP_data. + new_aquifer.time_constant_ = rst_aquifer.time_constant; + } +} Aquifetp Aquifetp::serializeObject() { diff --git a/src/opm/parser/eclipse/EclipseState/EclipseState.cpp b/src/opm/parser/eclipse/EclipseState/EclipseState.cpp index 2a9ad60f8..aa7279bee 100644 --- a/src/opm/parser/eclipse/EclipseState/EclipseState.cpp +++ b/src/opm/parser/eclipse/EclipseState/EclipseState.cpp @@ -26,6 +26,8 @@ #include #include +#include + #include #include #include @@ -46,7 +48,6 @@ #include #include - namespace Opm { @@ -266,6 +267,11 @@ AquiferConfig load_aquifers(const Deck& deck, const TableManager& tables, NNC& i return this->tracer_config; } + void EclipseState::loadRestartAquifers(const RestartIO::RstAquifer& aquifers) { + if (aquifers.hasAnalyticAquifers()) + this->aquifer_config.loadFromRestart(aquifers, this->m_tables); + } + void EclipseState::initTransMult() { const auto& fp = this->field_props; if (fp.has_double("MULTX")) this->m_transMult.applyMULT(fp.get_global_double("MULTX") , FaceDir::XPlus);