Merge pull request #2809 from bska/ext-network-restart

Support Loading Extended Network Model From Restart

CMakeLists_files.cmake

@@ -298,6 +298,7 @@ if(ENABLE_ECL_OUTPUT)
           src/opm/io/eclipse/rst/connection.cpp
           src/opm/io/eclipse/rst/group.cpp
           src/opm/io/eclipse/rst/header.cpp
+          src/opm/io/eclipse/rst/network.cpp
           src/opm/io/eclipse/rst/udq.cpp
           src/opm/io/eclipse/rst/segment.cpp
           src/opm/io/eclipse/rst/state.cpp
@@ -910,6 +911,7 @@ if(ENABLE_ECL_OUTPUT)
         opm/io/eclipse/rst/connection.hpp
         opm/io/eclipse/rst/group.hpp
         opm/io/eclipse/rst/header.hpp
+        opm/io/eclipse/rst/network.hpp
         opm/io/eclipse/rst/segment.hpp
         opm/io/eclipse/rst/state.hpp
         opm/io/eclipse/rst/udq.hpp

opm/io/eclipse/rst/network.hpp

@@ -0,0 +1,95 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef RST_NETWORK_HPP
+#define RST_NETWORK_HPP
+
+#include <memory>
+#include <optional>
+#include <string>
+#include <vector>
+
+namespace Opm {
+    class UnitSystem;
+} // namespace Opm
+
+namespace Opm { namespace EclIO {
+    class RestartFileView;
+}} // namespace Opm::EclIO
+
+namespace Opm { namespace RestartIO {
+
+    class RstNetwork
+    {
+    public:
+        /// Single branch in extended network model.
+        struct Branch
+        {
+            /// Downtree node.  Index into 'nodes' array.
+            int down{-1};
+
+            /// Uptree node.  Index into 'nodes' array.
+            int up{-1};
+
+            /// One-based VFP table ID.
+            int vfp{-1};
+        };
+
+        /// Single node in extended network model.
+        struct Node
+        {
+            /// Name of network node.
+            std::string name{};
+
+            /// Fixed pressure for terminal node.  Nullopt if not terminal.
+            std::optional<double> terminal_pressure{};
+
+            /// Group whose rate target the choking mechanism attempts to
+            /// match.  Nullopt if this node does not act as a choke or if
+            /// choking is disabled.
+            std::optional<std::string> as_choke{};
+
+            /// Whether or not to include lift gas of subordinate wells as
+            /// part of the produced gas entering the network at this node.
+            bool add_lift_gas{false};
+        };
+
+        explicit RstNetwork(std::shared_ptr<EclIO::RestartFileView> rstView,
+                            const UnitSystem&                       usys);
+
+        bool isActive() const;
+
+        const std::vector<Branch>& branches() const
+        {
+            return this->branches_;
+        }
+
+        const std::vector<Node>& nodes() const
+        {
+            return this->nodes_;
+        }
+
+    private:
+        std::vector<Branch> branches_{};
+        std::vector<Node> nodes_{};
+    };
+
+}} // namespace Opm::RestartIO
+
+#endif // RST_NETWORK_HPP

opm/io/eclipse/rst/state.hpp

@@ -24,22 +24,22 @@
 #include <unordered_map>
 #include <vector>
 
-#include <opm/io/eclipse/rst/header.hpp>
 #include <opm/io/eclipse/rst/aquifer.hpp>
-#include <opm/io/eclipse/rst/group.hpp>
-#include <opm/io/eclipse/rst/well.hpp>
-#include <opm/io/eclipse/rst/udq.hpp>
 #include <opm/io/eclipse/rst/action.hpp>
+#include <opm/io/eclipse/rst/group.hpp>
+#include <opm/io/eclipse/rst/header.hpp>
+#include <opm/io/eclipse/rst/network.hpp>
+#include <opm/io/eclipse/rst/udq.hpp>
+#include <opm/io/eclipse/rst/well.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Tuning.hpp>
 
 #include <opm/parser/eclipse/Units/UnitSystem.hpp>
 
-#include <opm/parser/eclipse/EclipseState/Schedule/Tuning.hpp>
-#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
-
 namespace Opm {
     class EclipseGrid;
     class Parser;
-    class Parser;
 } // namespace Opm
 
 namespace Opm { namespace EclIO {
@@ -62,6 +62,7 @@ struct RstState {
     ::Opm::UnitSystem unit_system;
     RstHeader header;
     RstAquifer aquifers;
+    RstNetwork network;
     std::vector<RstWell> wells;
     std::vector<RstGroup> groups;
     std::vector<RstUDQ> udqs;

src/opm/io/eclipse/rst/network.cpp

@@ -0,0 +1,250 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/io/eclipse/rst/network.hpp>
+
+#include <opm/io/eclipse/RestartFileView.hpp>
+
+#include <opm/output/eclipse/VectorItems/intehead.hpp>
+#include <opm/output/eclipse/VectorItems/network.hpp>
+
+#include <opm/parser/eclipse/Units/UnitSystem.hpp>
+
+#include <cstddef>
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <boost/range.hpp>
+
+namespace VI = ::Opm::RestartIO::Helpers::VectorItems;
+
+namespace {
+    template <typename T>
+    boost::iterator_range<typename std::vector<T>::const_iterator>
+    getDataWindow(const std::vector<T>& arr,
+                  const std::size_t     windowSize,
+                  const std::size_t     entity)
+    {
+        const auto off = windowSize * entity;
+
+        auto begin = arr.begin() + off;
+        auto end   = begin       + windowSize;
+
+        return { begin, end };
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+class BranchVectors
+{
+public:
+    template <typename T>
+    using Window = boost::iterator_range<
+        typename std::vector<T>::const_iterator
+    >;
+
+    explicit BranchVectors(const std::vector<int>&                      intehead,
+                           std::shared_ptr<Opm::EclIO::RestartFileView> rst_view);
+
+    std::size_t numActiveBranches() const
+    {
+        return this->numActiveBranches_;
+    }
+
+    Window<int> ibran(const std::size_t branchID) const;
+
+private:
+    std::size_t numActiveBranches_;
+    std::size_t numIBranElem_;
+
+    std::shared_ptr<Opm::EclIO::RestartFileView> rstView_;
+};
+
+BranchVectors::BranchVectors(const std::vector<int>&                      intehead,
+                             std::shared_ptr<Opm::EclIO::RestartFileView> rst_view)
+    : numActiveBranches_(intehead[VI::intehead::NOACTBR])
+    , numIBranElem_     (intehead[VI::intehead::NIBRAN])
+    , rstView_          (std::move(rst_view))
+{}
+
+BranchVectors::Window<int>
+BranchVectors::ibran(const std::size_t branchID) const
+{
+    return getDataWindow(this->rstView_->getKeyword<int>("IBRAN"),
+                         this->numIBranElem_, branchID);
+}
+
+// ---------------------------------------------------------------------------
+
+class NodeVectors
+{
+public:
+    template <typename T>
+    using Window = boost::iterator_range<
+        typename std::vector<T>::const_iterator
+    >;
+
+    explicit NodeVectors(const std::vector<int>&                      intehead,
+                         std::shared_ptr<Opm::EclIO::RestartFileView> rst_view);
+
+    std::size_t numActiveNodes() const
+    {
+        return this->numActiveNodes_;
+    }
+
+    Window<int>         inode(const std::size_t nodeID) const;
+    Window<double>      rnode(const std::size_t nodeID) const;
+    Window<std::string> znode(const std::size_t nodeID) const;
+
+private:
+    std::size_t numActiveNodes_;
+    std::size_t numINodeElem_;
+    std::size_t numRNodeElem_;
+    std::size_t numZNodeElem_;
+
+    std::shared_ptr<Opm::EclIO::RestartFileView> rstView_;
+};
+
+NodeVectors::NodeVectors(const std::vector<int>&                      intehead,
+                         std::shared_ptr<Opm::EclIO::RestartFileView> rst_view)
+    : numActiveNodes_(intehead[VI::intehead::NOACTNOD])
+    , numINodeElem_  (intehead[VI::intehead::NINODE])
+    , numRNodeElem_  (intehead[VI::intehead::NRNODE])
+    , numZNodeElem_  (intehead[VI::intehead::NZNODE])
+    , rstView_       (std::move(rst_view))
+{}
+
+NodeVectors::Window<int>
+NodeVectors::inode(const std::size_t nodeID) const
+{
+    return getDataWindow(this->rstView_->getKeyword<int>("INODE"),
+                         this->numINodeElem_, nodeID);
+}
+
+NodeVectors::Window<double>
+NodeVectors::rnode(const std::size_t nodeID) const
+{
+    return getDataWindow(this->rstView_->getKeyword<double>("RNODE"),
+                         this->numRNodeElem_, nodeID);
+}
+
+NodeVectors::Window<std::string>
+NodeVectors::znode(const std::size_t nodeID) const
+{
+    return getDataWindow(this->rstView_->getKeyword<std::string>("ZNODE"),
+                         this->numZNodeElem_, nodeID);
+}
+
+// ---------------------------------------------------------------------------
+
+namespace {
+    int branch_vfp_no_pressure_loss()
+    {
+        return 9999;
+    }
+
+    template <typename IBranArray>
+    Opm::RestartIO::RstNetwork::Branch make_branch(const IBranArray ibran)
+    {
+        using Ix = VI::IBran::index;
+
+        auto branch = Opm::RestartIO::RstNetwork::Branch{};
+
+        branch.down = ibran[Ix::DownTreeNode] - 1; // Index into ZNODE
+        branch.up   = ibran[Ix::UpTreeNode]   - 1; // Index into ZNODE
+        branch.vfp  = (ibran[Ix::VfpTableNo] > 0)
+            ? ibran[Ix::VfpTableNo] // One-based VFP table ID
+            : branch_vfp_no_pressure_loss();
+
+        return branch;
+    }
+
+    template <typename INodeArray, typename RNodeArray>
+    Opm::RestartIO::RstNetwork::Node
+    make_node(const std::string&     name,
+              const Opm::UnitSystem& usys,
+              const INodeArray       inode,
+              const RNodeArray       rnode)
+    {
+        auto node = Opm::RestartIO::RstNetwork::Node{};
+        node.name = name;
+
+        if (inode[VI::INode::index::FixedPresNode] != 0) {
+            node.terminal_pressure =
+                usys.to_si(Opm::UnitSystem::measure::pressure,
+                           rnode[VI::RNode::index::PressureLimit]);
+        }
+
+        return node;
+    }
+
+    std::vector<Opm::RestartIO::RstNetwork::Branch>
+    create_branches(std::shared_ptr<Opm::EclIO::RestartFileView> rstView)
+    {
+        auto branches = std::vector<Opm::RestartIO::RstNetwork::Branch>{};
+
+        const auto branchData = BranchVectors { rstView->intehead(), rstView };
+        const auto numBranches = branchData.numActiveBranches();
+
+        branches.reserve(numBranches);
+
+        for (auto branchID = 0*numBranches; branchID < numBranches; ++branchID) {
+            branches.push_back(make_branch(branchData.ibran(branchID)));
+        }
+
+        return branches;
+    }
+
+    std::vector<Opm::RestartIO::RstNetwork::Node>
+    create_nodes(std::shared_ptr<Opm::EclIO::RestartFileView> rstView,
+                 const Opm::UnitSystem&                       usys)
+    {
+        auto nodes = std::vector<Opm::RestartIO::RstNetwork::Node>{};
+
+        const auto nodeData = NodeVectors { rstView->intehead(), rstView };
+        const auto numNodes = nodeData.numActiveNodes();
+
+        nodes.reserve(numNodes);
+
+        for (auto nodeID = 0*numNodes; nodeID < numNodes; ++nodeID) {
+            const auto znode = nodeData.znode(nodeID);
+
+            nodes.push_back(make_node(znode[VI::ZNode::index::NodeName], usys,
+                                      nodeData.inode(nodeID),
+                                      nodeData.rnode(nodeID)));
+        }
+
+        return nodes;
+    }
+}
+
+Opm::RestartIO::RstNetwork::RstNetwork(std::shared_ptr<EclIO::RestartFileView> rstView,
+                                       const UnitSystem&                       usys)
+    : branches_{ create_branches(rstView) }
+    , nodes_   { create_nodes   (rstView, usys) }
+{}
+
+bool Opm::RestartIO::RstNetwork::isActive() const
+{
+    return ! (this->branches_.empty() ||
+              this->nodes_.empty());
+}

src/opm/io/eclipse/rst/state.cpp

@@ -16,33 +16,41 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <opm/io/eclipse/rst/state.hpp>
+
+#include <opm/io/eclipse/RestartFileView.hpp>
+
+#include <opm/io/eclipse/rst/aquifer.hpp>
+#include <opm/io/eclipse/rst/action.hpp>
+#include <opm/io/eclipse/rst/connection.hpp>
+#include <opm/io/eclipse/rst/group.hpp>
+#include <opm/io/eclipse/rst/header.hpp>
+#include <opm/io/eclipse/rst/network.hpp>
+#include <opm/io/eclipse/rst/udq.hpp>
+#include <opm/io/eclipse/rst/well.hpp>
+
+#include <opm/common/utility/String.hpp>
+#include <opm/common/utility/TimeService.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Schedule/Action/Actdims.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Action/Condition.hpp>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+
+#include <opm/output/eclipse/UDQDims.hpp>
+#include <opm/output/eclipse/VectorItems/connection.hpp>
+#include <opm/output/eclipse/VectorItems/doubhead.hpp>
+#include <opm/output/eclipse/VectorItems/intehead.hpp>
+#include <opm/output/eclipse/VectorItems/well.hpp>
+#include <opm/output/eclipse/WriteRestartHelpers.hpp>
+
 #include <algorithm>
 #include <iterator>
 #include <memory>
 #include <numeric>
 #include <optional>
 
-#include <opm/io/eclipse/RestartFileView.hpp>
-
-#include <opm/io/eclipse/rst/header.hpp>
-#include <opm/io/eclipse/rst/connection.hpp>
-#include <opm/io/eclipse/rst/well.hpp>
-#include <opm/io/eclipse/rst/state.hpp>
-#include <opm/output/eclipse/WriteRestartHelpers.hpp>
-
-#include <opm/common/utility/TimeService.hpp>
-#include <opm/parser/eclipse/Deck/Deck.hpp>
-#include <opm/common/utility/String.hpp>
-#include <opm/parser/eclipse/Parser/Parser.hpp>
-#include <opm/output/eclipse/UDQDims.hpp>
-#include <opm/output/eclipse/VectorItems/connection.hpp>
-#include <opm/output/eclipse/VectorItems/well.hpp>
-#include <opm/output/eclipse/VectorItems/intehead.hpp>
-#include <opm/output/eclipse/VectorItems/doubhead.hpp>
-#include <opm/parser/eclipse/EclipseState/Schedule/Action/Actdims.hpp>
-#include <opm/parser/eclipse/EclipseState/Schedule/Action/Condition.hpp>
-
-
 namespace {
     std::string
     udq_define(const std::vector<std::string>& zudl,
@@ -78,6 +86,7 @@ RstState::RstState(std::shared_ptr<EclIO::RestartFileView> rstView,
     : unit_system(rstView->intehead()[VI::intehead::UNIT])
     , header(unit_system, rstView->intehead(), rstView->logihead(), rstView->doubhead())
     , aquifers(rstView, grid, unit_system)
+    , network(rstView, unit_system)
 {
     this->load_tuning(rstView->intehead(), rstView->doubhead());
 }

src/opm/parser/eclipse/EclipseState/Schedule/Schedule.cpp

@@ -47,6 +47,7 @@
 #include <opm/parser/eclipse/Deck/DeckSection.hpp>
 #include <opm/parser/eclipse/Parser/ErrorGuard.hpp>
 #include <opm/parser/eclipse/Parser/ParseContext.hpp>
+#include <opm/parser/eclipse/Parser/ParserKeywords/B.hpp>
 #include <opm/parser/eclipse/Parser/ParserKeywords/E.hpp>
 #include <opm/parser/eclipse/Parser/ParserKeywords/P.hpp>
 #include <opm/parser/eclipse/Parser/ParserKeywords/R.hpp>
@@ -1561,6 +1562,45 @@ namespace {
 
         if (!rst_state.wlists.empty())
             this->snapshots.back().wlist_manager.update( WListManager(rst_state) );
+
+        if (rst_state.network.isActive()) {
+            auto network = this->snapshots.back().network();
+
+            // Note: We presently support only the default value of BRANPROP(4).
+            const auto alq_value =
+                ParserKeywords::BRANPROP::ALQ::defaultValue;
+
+            const auto& rst_nodes = rst_state.network.nodes();
+            for (const auto& rst_branch : rst_state.network.branches()) {
+                if ((rst_branch.down < 0) || (rst_branch.up < 0)) {
+                    // Prune branches to non-existent nodes.
+                    continue;
+                }
+
+                const auto& downtree_node = rst_nodes[rst_branch.down].name;
+                const auto& uptree_node = rst_nodes[rst_branch.up].name;
+
+                network.add_branch({ downtree_node, uptree_node, rst_branch.vfp, alq_value });
+            }
+
+            for (const auto& rst_node : rst_nodes) {
+                auto node = Network::Node { rst_node.name };
+
+                if (rst_node.terminal_pressure.has_value()) {
+                    node.terminal_pressure(rst_node.terminal_pressure.value());
+                }
+
+                if (rst_node.as_choke.has_value()) {
+                    node.as_choke(rst_node.as_choke.value());
+                }
+
+                node.add_gas_lift_gas(rst_node.add_lift_gas);
+
+                network.add_node(std::move(node));
+            }
+
+            this->snapshots.back().network.update(std::move(network));
+        }
     }
 
     std::shared_ptr<const Python> Schedule::python() const

tests/parser/ScheduleRestartTests.cpp

@@ -123,7 +123,7 @@ std::tuple<Schedule, Schedule, RestartIO::RstState> load_schedule_pair(const std
     EclipseState ecl_state_restart(restart_deck);
     Schedule restart_sched(restart_deck, ecl_state_restart, python, {}, &rst_state);
 
-    return {sched, restart_sched, rst_state};
+    return {sched, restart_sched, std::move(rst_state)};
 }