diff --git a/opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp index b7304505a6..74ebdd64cc 100644 --- a/opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp @@ -29,7 +29,6 @@ #include #include -#include #include #include @@ -144,10 +143,10 @@ updatePressures(const int reportStepIdx, const auto previous_node_pressures = node_pressures_; - node_pressures_ = well_model_.groupStateHelper(). - computeNetworkPressures(network, - *well_model_.getVFPProperties().getProd(), - well_model_.comm()); + node_pressures_ = this->computePressures(network, + *well_model_.getVFPProperties().getProd(), + reportStepIdx, + well_model_.comm()); // here, the network imbalance is the difference between the previous nodal pressure and the new nodal pressure Scalar network_imbalance = 0.; @@ -233,11 +232,10 @@ assignNodeValues(std::map& nodevalues, return; } - auto converged_pressures = - well_model_.groupStateHelper(). - computeNetworkPressures(network, - *well_model_.getVFPProperties().getProd(), - well_model_.comm()); + auto converged_pressures = this->computePressures(network, + *well_model_.getVFPProperties().getProd(), + reportStepIdx, + well_model_.comm()); for (const auto& [node, converged_pressure] : converged_pressures) { auto it = nodevalues.find(node); assert(it != nodevalues.end() ); @@ -285,6 +283,169 @@ initializeWell(WellInterfaceGeneric& well) } } +template +std::map +BlackoilWellModelNetworkGeneric:: +computePressures(const Network::ExtNetwork& network, + const VFPProdProperties& vfp_prod_props, + const int reportStepIdx, + const Parallel::Communication& comm) const +{ + // TODO: Only dealing with production networks for now. + OPM_TIMEFUNCTION(); + if (!network.active()) { + return {}; + } + + const auto& group_state = well_model_.groupStateHelper().groupState(); + const auto& well_state = well_model_.groupStateHelper().wellState(); + + std::map node_pressures; + auto roots = network.roots(); + for (const auto& root : roots) { + // Fixed pressure nodes of the network are the roots of trees. + // Leaf nodes must correspond to groups in the group structure. + // Let us first find all leaf nodes of the network. We also + // create a vector of all nodes, ordered so that a child is + // always after its parent. + std::stack children; + std::set leaf_nodes; + std::vector root_to_child_nodes; + // children.push(network.root().name()); + children.push(root.get().name()); + while (!children.empty()) { + const auto node = children.top(); + children.pop(); + root_to_child_nodes.push_back(node); + auto branches = network.downtree_branches(node); + if (branches.empty()) { + leaf_nodes.insert(node); + } + for (const auto& branch : branches) { + children.push(branch.downtree_node()); + } + } + assert(children.empty()); + + // Starting with the leaf nodes of the network, get the flow rates + // from the corresponding groups. + std::map> node_inflows; + const std::vector zero_rates(3, 0.0); + for (const auto& node : leaf_nodes) { + // Guard against empty leaf nodes (may not be present in GRUPTREE) + if (!well_model_.groupStateHelper().groupState().has_production_rates(node)) { + node_inflows[node] = zero_rates; + continue; + } + + node_inflows[node] = group_state.network_leaf_node_production_rates(node); + // Add the ALQ amounts to the gas rates if requested. + if (network.node(node).add_gas_lift_gas()) { + const auto& group = well_model_.schedule().getGroup(node, reportStepIdx); + Scalar alq = 0.0; + for (const std::string& wellname : group.wells()) { + const Well& well = well_model_.schedule().getWell(wellname, reportStepIdx); + if (well.isInjector() || !well_state.isOpen(wellname)) + continue; + const Scalar efficiency = well.getEfficiencyFactor(/*network*/ true) + * well_state.getGlobalEfficiencyScalingFactor(wellname); + const auto& well_index = well_state.index(wellname); + if (well_index.has_value() && + well_state.wellIsOwned(well_index.value(), wellname)) + { + alq += well_state.well(wellname).alq_state.get() * efficiency; + } + } + alq = comm.sum(alq); + node_inflows[node][IndexTraits::gasPhaseIdx] += alq; + } + } + + // Accumulate in the network, towards the roots. Note that a + // root (i.e. fixed pressure node) can still be contributing + // flow towards other nodes in the network, i.e. a node is + // the root of a subtree. + auto child_to_root_nodes = root_to_child_nodes; + std::reverse(child_to_root_nodes.begin(), child_to_root_nodes.end()); + for (const auto& node : child_to_root_nodes) { + const auto upbranch = network.uptree_branch(node); + if (upbranch) { + // Add downbranch rates to upbranch. + std::vector& up = node_inflows[(*upbranch).uptree_node()]; + const std::vector& down = node_inflows[node]; + // We now also support NEFAC + const Scalar efficiency = network.node(node).efficiency(); + if (up.empty()) { + up = std::vector(down.size(), 0.0); + } + assert(up.size() == down.size()); + for (std::size_t ii = 0; ii < up.size(); ++ii) { + up[ii] += efficiency * down[ii]; + } + } + } + + // Going the other way (from roots to leafs), calculate the pressure + // at each node using VFP tables and rates. + // std::map node_pressures; + for (const auto& node : root_to_child_nodes) { + auto press = network.node(node).terminal_pressure(); + if (press) { + node_pressures[node] = *press; + } else { + const auto upbranch = network.uptree_branch(node); + assert(upbranch); + const Scalar up_press = node_pressures[(*upbranch).uptree_node()]; + const auto vfp_table = (*upbranch).vfp_table(); + if (vfp_table) { + OPM_TIMEBLOCK(NetworkVfpCalculations); + // The rates are here positive, but the VFP code expects the + // convention that production rates are negative, so we must + // take a copy and flip signs. + auto rates = node_inflows[node]; + std::transform( + rates.begin(), rates.end(), rates.begin(), [](const auto r) { return -r; }); + assert(rates.size() == 3); + // NB! ALQ in extended network is never implicitly the gas lift rate (GRAT), i.e., the + // gas lift rates only enters the network pressure calculations through the rates + // (e.g., in GOR calculations) unless a branch ALQ is set in BRANPROP. + // + // @TODO: Standard network + Scalar alq = (*upbranch).alq_value().value_or(0.0); + node_pressures[node] + = vfp_prod_props.bhp(*vfp_table, + rates[IndexTraits::waterPhaseIdx], + rates[IndexTraits::oilPhaseIdx], + rates[IndexTraits::gasPhaseIdx], + up_press, + alq, + 0.0, // explicit_wfr + 0.0, // explicit_gfr + false); // use_expvfp we dont support explicit lookup +#define EXTRA_DEBUG_NETWORK 0 +#if EXTRA_DEBUG_NETWORK + std::ostringstream oss; + oss << "parent: " << (*upbranch).uptree_node() << " child: " << node << " rates = [ " + << rates[0] * 86400 << ", " << rates[1] * 86400 << ", " << rates[2] * 86400 << " ]" + << " p(parent) = " << up_press / 1e5 << " p(child) = " << node_pressures[node] / 1e5 + << std::endl; + OpmLog::debug(oss.str()); +#endif + } else { + // Table number specified as 9999 in the deck, no pressure loss. + if (network.node(node).as_choke()) { + // Node pressure is set to the group THP. + node_pressures[node] = well_model_.groupStateHelper().groupState().well_group_thp(node); + } else { + node_pressures[node] = up_press; + } + } + } + } + } + return node_pressures; +} + template bool BlackoilWellModelNetworkGeneric:: operator==(const BlackoilWellModelNetworkGeneric& rhs) const diff --git a/opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp index 56083315f0..95989a6d25 100644 --- a/opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp @@ -23,7 +23,10 @@ #ifndef OPM_BLACKOILWELLMODEL_NETWORK_GENERIC_HEADER_INCLUDED #define OPM_BLACKOILWELLMODEL_NETWORK_GENERIC_HEADER_INCLUDED +#include + #include + #include #include @@ -34,6 +37,7 @@ namespace Opm { class Schedule; template class BlackoilWellModelGeneric; template class WellInterfaceGeneric; + template class VFPProdProperties; } namespace Opm { @@ -100,6 +104,13 @@ public: bool operator==(const BlackoilWellModelNetworkGeneric& rhs) const; protected: + std::map + computePressures(const Network::ExtNetwork& network, + const VFPProdProperties& vfp_prod_props, + const int reportStepIdx, + const Parallel::Communication& comm) const; + + bool active_{false}; BlackoilWellModelGeneric& well_model_; diff --git a/opm/simulators/wells/GroupStateHelper.cpp b/opm/simulators/wells/GroupStateHelper.cpp index 93f0b77c4f..cddb2cc15f 100644 --- a/opm/simulators/wells/GroupStateHelper.cpp +++ b/opm/simulators/wells/GroupStateHelper.cpp @@ -419,163 +419,6 @@ GroupStateHelper::checkGroupConstraintsProd(const std::stri return std::make_pair(current_rate_available > target_rate_available, scale); } -template -std::map -GroupStateHelper::computeNetworkPressures(const Network::ExtNetwork& network, - const VFPProdProperties& vfp_prod_props, - const Parallel::Communication& comm) const -{ - // TODO: Only dealing with production networks for now. - OPM_TIMEFUNCTION(); - if (!network.active()) { - return {}; - } - - std::map node_pressures; - auto roots = network.roots(); - for (const auto& root : roots) { - // Fixed pressure nodes of the network are the roots of trees. - // Leaf nodes must correspond to groups in the group structure. - // Let us first find all leaf nodes of the network. We also - // create a vector of all nodes, ordered so that a child is - // always after its parent. - std::stack children; - std::set leaf_nodes; - std::vector root_to_child_nodes; - // children.push(network.root().name()); - children.push(root.get().name()); - while (!children.empty()) { - const auto node = children.top(); - children.pop(); - root_to_child_nodes.push_back(node); - auto branches = network.downtree_branches(node); - if (branches.empty()) { - leaf_nodes.insert(node); - } - for (const auto& branch : branches) { - children.push(branch.downtree_node()); - } - } - assert(children.empty()); - - // Starting with the leaf nodes of the network, get the flow rates - // from the corresponding groups. - std::map> node_inflows; - const std::vector zero_rates(3, 0.0); - for (const auto& node : leaf_nodes) { - // Guard against empty leaf nodes (may not be present in GRUPTREE) - if (!this->groupState().has_production_rates(node)) { - node_inflows[node] = zero_rates; - continue; - } - - node_inflows[node] = this->groupState().network_leaf_node_production_rates(node); - // Add the ALQ amounts to the gas rates if requested. - if (network.node(node).add_gas_lift_gas()) { - const auto& group = this->schedule_.getGroup(node, this->report_step_); - Scalar alq = 0.0; - for (const std::string& wellname : group.wells()) { - const Well& well = this->schedule_.getWell(wellname, this->report_step_); - if (well.isInjector() || !this->wellState().isOpen(wellname)) - continue; - const Scalar efficiency = well.getEfficiencyFactor(/*network*/ true) - * this->wellState().getGlobalEfficiencyScalingFactor(wellname); - const auto& well_index = this->wellState().index(wellname); - if (well_index.has_value() - && this->wellState().wellIsOwned(well_index.value(), wellname)) { - alq += this->wellState().well(wellname).alq_state.get() * efficiency; - } - } - alq = comm.sum(alq); - node_inflows[node][IndexTraits::gasPhaseIdx] += alq; - } - } - - // Accumulate in the network, towards the roots. Note that a - // root (i.e. fixed pressure node) can still be contributing - // flow towards other nodes in the network, i.e. a node is - // the root of a subtree. - auto child_to_root_nodes = root_to_child_nodes; - std::reverse(child_to_root_nodes.begin(), child_to_root_nodes.end()); - for (const auto& node : child_to_root_nodes) { - const auto upbranch = network.uptree_branch(node); - if (upbranch) { - // Add downbranch rates to upbranch. - std::vector& up = node_inflows[(*upbranch).uptree_node()]; - const std::vector& down = node_inflows[node]; - // We now also support NEFAC - const Scalar efficiency = network.node(node).efficiency(); - if (up.empty()) { - up = std::vector(down.size(), 0.0); - } - assert(up.size() == down.size()); - for (std::size_t ii = 0; ii < up.size(); ++ii) { - up[ii] += efficiency * down[ii]; - } - } - } - - // Going the other way (from roots to leafs), calculate the pressure - // at each node using VFP tables and rates. - // std::map node_pressures; - for (const auto& node : root_to_child_nodes) { - auto press = network.node(node).terminal_pressure(); - if (press) { - node_pressures[node] = *press; - } else { - const auto upbranch = network.uptree_branch(node); - assert(upbranch); - const Scalar up_press = node_pressures[(*upbranch).uptree_node()]; - const auto vfp_table = (*upbranch).vfp_table(); - if (vfp_table) { - OPM_TIMEBLOCK(NetworkVfpCalculations); - // The rates are here positive, but the VFP code expects the - // convention that production rates are negative, so we must - // take a copy and flip signs. - auto rates = node_inflows[node]; - std::transform( - rates.begin(), rates.end(), rates.begin(), [](const auto r) { return -r; }); - assert(rates.size() == 3); - // NB! ALQ in extended network is never implicitly the gas lift rate (GRAT), i.e., the - // gas lift rates only enters the network pressure calculations through the rates - // (e.g., in GOR calculations) unless a branch ALQ is set in BRANPROP. - // - // @TODO: Standard network - Scalar alq = (*upbranch).alq_value().value_or(0.0); - node_pressures[node] - = vfp_prod_props.bhp(*vfp_table, - rates[IndexTraits::waterPhaseIdx], - rates[IndexTraits::oilPhaseIdx], - rates[IndexTraits::gasPhaseIdx], - up_press, - alq, - 0.0, // explicit_wfr - 0.0, // explicit_gfr - false); // use_expvfp we dont support explicit lookup -#define EXTRA_DEBUG_NETWORK 0 -#if EXTRA_DEBUG_NETWORK - std::ostringstream oss; - oss << "parent: " << (*upbranch).uptree_node() << " child: " << node << " rates = [ " - << rates[0] * 86400 << ", " << rates[1] * 86400 << ", " << rates[2] * 86400 << " ]" - << " p(parent) = " << up_press / 1e5 << " p(child) = " << node_pressures[node] / 1e5 - << std::endl; - OpmLog::debug(oss.str()); -#endif - } else { - // Table number specified as 9999 in the deck, no pressure loss. - if (network.node(node).as_choke()) { - // Node pressure is set to the group THP. - node_pressures[node] = this->groupState().well_group_thp(node); - } else { - node_pressures[node] = up_press; - } - } - } - } - } - return node_pressures; -} - template Scalar GroupStateHelper::getGuideRate(const std::string& name, diff --git a/opm/simulators/wells/GroupStateHelper.hpp b/opm/simulators/wells/GroupStateHelper.hpp index 6348c50988..d9b1664588 100644 --- a/opm/simulators/wells/GroupStateHelper.hpp +++ b/opm/simulators/wells/GroupStateHelper.hpp @@ -144,10 +144,6 @@ public: const bool check_guide_rate, DeferredLogger& deferred_logger) const; - std::map computeNetworkPressures(const Network::ExtNetwork& network, - const VFPProdProperties& vfp_prod_props, - const Parallel::Communication& comm) const; - Scalar getGuideRate(const std::string& name, const GuideRateModel::Target target) const; GuideRate::RateVector getProductionGroupRateVector(const std::string& group_name) const;