diff --git a/opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp index b461391e2c8..ec366070bc7 100644 --- a/opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp @@ -35,10 +35,70 @@ #include #include +#include + +#include + +#include +#include + +#include #include namespace Opm { +namespace details { + /// Helper to check if any network (production, gas injection, water injection) is active at a given time step. + bool anyNetworkActive(const Schedule& schedule, const int timeStepIdx) + { + const auto& sstate = schedule[timeStepIdx]; + return sstate.network().active() + || (sstate.injectionNetwork.get_ptr(Phase::GAS) != nullptr + && sstate.injectionNetwork.get_ptr(Phase::GAS)->active()) + || (sstate.injectionNetwork.get_ptr(Phase::WATER) != nullptr + && sstate.injectionNetwork.get_ptr(Phase::WATER)->active()); + } + + /// Helper to get all active networks (production, gas injection, water injection) at a given time step. + std::vector + activeNetworks(const Schedule& schedule, const int timeStepIdx) + { + std::vector active_networks; + const auto& sstate = schedule[timeStepIdx]; + if (sstate.network().active()) { + active_networks.push_back({NetworkDomain::Production, std::cref(sstate.network())}); + } + if (sstate.injectionNetwork.get_ptr(Phase::GAS) != nullptr + && sstate.injectionNetwork.get_ptr(Phase::GAS)->active()) { + active_networks.push_back({NetworkDomain::InjectionGas, std::cref(*sstate.injectionNetwork.get_ptr(Phase::GAS))}); + } + if (sstate.injectionNetwork.get_ptr(Phase::WATER) != nullptr + && sstate.injectionNetwork.get_ptr(Phase::WATER)->active()) { + active_networks.push_back({NetworkDomain::InjectionWater, std::cref(*sstate.injectionNetwork.get_ptr(Phase::WATER))}); + } + return active_networks; + } + + /// Helper to get all networks (production, gas injection, water injection) at a given time step, + /// whether active or not. + std::vector + networks(const Schedule& schedule, const int timeStepIdx) + { + std::vector nw; + const auto& sstate = schedule[timeStepIdx]; + nw.push_back({NetworkDomain::Production, std::cref(sstate.network())}); + if (sstate.injectionNetwork.get_ptr(Phase::GAS) != nullptr) { + nw.push_back({NetworkDomain::InjectionGas, std::cref(*sstate.injectionNetwork.get_ptr(Phase::GAS))}); + } + if (sstate.injectionNetwork.get_ptr(Phase::WATER) != nullptr) { + nw.push_back({NetworkDomain::InjectionWater, std::cref(*sstate.injectionNetwork.get_ptr(Phase::WATER))}); + } + return nw; + } + +} // namespace details + + template BlackoilWellModelNetworkGeneric:: BlackoilWellModelNetworkGeneric(BlackoilWellModelGeneric& well_model) @@ -59,6 +119,7 @@ setFromRestart(const std::optional>& node_pressure this->node_pressures_[it.first] = it.second; } } + this->syncProductionDomainState_(); } } @@ -66,12 +127,22 @@ template void BlackoilWellModelNetworkGeneric:: updateActiveState(const int report_step) { - const auto& network = well_model_.schedule()[report_step].network(); + this->active_ = false; + for (const auto& network : details::activeNetworks(well_model_.schedule(), report_step)) { + updateActiveStateImpl(network.network.get()); + } + this->active_ = well_model_.comm().max(active_); +} + +template +void BlackoilWellModelNetworkGeneric:: +updateActiveStateImpl(const Network::ExtNetwork& network) +{ + // Accumulates into active_ across the domains; an inactive network must not + // clear what an earlier domain set. if (!network.active()) { - this->active_ = false; return; } - bool network_active = false; for (const auto& well : well_model_.genericWells()) { const bool is_partof_network = network.has_node(well->wellEcl().groupName()); @@ -108,17 +179,22 @@ updateActiveState(const int report_step) } } #endif - this->active_ = well_model_.comm().max(network_active); + this->active_ = this->active_ || network_active; } template bool BlackoilWellModelNetworkGeneric:: needPreStepRebalance(const int report_step) const { - const auto& network = well_model_.schedule()[report_step].network(); + const auto active_networks = details::activeNetworks(well_model_.schedule(), report_step); bool network_rebalance_necessary = false; for (const auto& well : well_model_.genericWells()) { - const bool is_partof_network = network.has_node(well->wellEcl().groupName()); + const bool is_partof_network = std::any_of(active_networks.begin(), + active_networks.end(), + [&](const auto& network) + { + return network.network.get().has_node(well->wellEcl().groupName()); + }); // TODO: we might find more relevant events to be included here (including network change events?) const auto& events = well_model_.wellState().well(well->indexOfWell()).events; if (is_partof_network && events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE)) { @@ -135,8 +211,7 @@ bool BlackoilWellModelNetworkGeneric:: shouldBalance(const int reportStepIdx) const { // if network is not active, we do not need to balance the network - const auto& network = well_model_.schedule()[reportStepIdx].network(); - if (!network.active()) { + if (!details::anyNetworkActive(well_model_.schedule(), reportStepIdx)) { return false; } @@ -161,11 +236,11 @@ bool BlackoilWellModelNetworkGeneric:: willBalanceOnNextIteration(const int reportStepIdx) const { // if network is not active, we do not need to balance the network - const auto& schedule_state = well_model_.schedule()[reportStepIdx]; - if (!schedule_state.network().active()) { + if (!details::anyNetworkActive(well_model_.schedule(), reportStepIdx)) { return false; } + const auto& schedule_state = well_model_.schedule()[reportStepIdx]; if (schedule_state.network_balance().mode() == Network::Balance::CalcMode::NUPCOL) { const int nupcol = schedule_state.nupcol(); return well_model_.iterationContext().withinNupcol(nupcol - 1); // Note the -1 here! @@ -184,19 +259,38 @@ updatePressures(const int reportStepIdx, const Scalar upper_update_bound) { OPM_TIMEFUNCTION(); - // Get the network and return if inactive (no wells in network at this time) - const auto& network = well_model_.schedule()[reportStepIdx].network(); - if (!network.active()) { - return 0.0; + for (auto& invalid_nodes : this->domain_invalid_nodes_) { + invalid_nodes.clear(); } + this->invalid_nodes_report_step_ = -1; - const auto previous_node_pressures = node_pressures_; + if (!details::anyNetworkActive(well_model_.schedule(), reportStepIdx)) { + return 0.0; + } - std::tie(node_pressures_, branch_data_) = this->computePressures(network, - *well_model_.getVFPProperties().getProd(), - well_model_.schedule().getUnits(), - reportStepIdx, - well_model_.comm()); + this->syncProductionDomainState_(); + const auto previous_node_pressures = this->domain_node_pressures_; + + for (const auto& network : details::activeNetworks(well_model_.schedule(), reportStepIdx)) { + NetworkPressures result; + if (network.domain == details::NetworkDomain::Production) { + result = this->computePressures(network.network.get(), + *well_model_.getVFPProperties().getProd(), + well_model_.schedule().getUnits(), + reportStepIdx, + well_model_.comm()); + } else { + result = this->computePressures(network.network.get(), + *well_model_.getVFPProperties().getInj(), + well_model_.schedule().getUnits(), + reportStepIdx, + well_model_.comm()); + } + this->nodePressures(network.domain) = std::move(result.node_pressures); + this->branchData(network.domain) = std::move(result.branch_data); + this->invalidNodes(network.domain) = std::move(result.invalid_nodes); + } + this->syncLegacyProductionState_(); // here, the network imbalance is the difference between the previous nodal pressure and the new nodal pressure Scalar network_imbalance = 0.; @@ -204,54 +298,114 @@ updatePressures(const int reportStepIdx, return network_imbalance; } - if (!previous_node_pressures.empty()) { - for (const auto& [name, new_pressure]: node_pressures_) { - if (previous_node_pressures.count(name) <= 0) { - if (std::abs(new_pressure) > network_imbalance) { - network_imbalance = std::abs(new_pressure); - } - continue; + for (const auto& network : details::activeNetworks(well_model_.schedule(), reportStepIdx)) { + auto& domain_pressures = this->nodePressures(network.domain); + const auto& invalid = this->invalidNodes(network.domain); + const auto& previous_domain_pressures = previous_node_pressures[details::domainIndex(network.domain)]; + + if (!invalid.empty()) { + // The VFP tables gave no pressure for these nodes (rate/pressure outside what + // the tables can deliver). Keep the previous value and report the network as + // unbalanced so that the wells get another chance to move into range. + network_imbalance = std::max(network_imbalance, upper_update_bound); + if (this->invalid_nodes_report_step_ != reportStepIdx) { + this->invalid_nodes_report_step_ = reportStepIdx; + OpmLog::warning(fmt::format("Network: no VFP solution for node(s) {} at report step {}; " + "keeping the previous node pressure(s).", + fmt::join(invalid, ", "), reportStepIdx + 1)); } - const auto pressure = previous_node_pressures.at(name); - const Scalar change = (new_pressure - pressure); - if (std::abs(change) > network_imbalance) { - network_imbalance = std::abs(change); + } + + if (!previous_domain_pressures.empty()) { + for (auto& [name, new_pressure]: domain_pressures) { + if (previous_domain_pressures.count(name) <= 0) { + if (std::abs(new_pressure) > network_imbalance) { + network_imbalance = std::abs(new_pressure); + } + continue; + } + + const auto pressure = previous_domain_pressures.at(name); + if (invalid.count(name) > 0) { + new_pressure = pressure; + continue; + } + const Scalar change = (new_pressure - pressure); + if (std::abs(change) > network_imbalance) { + network_imbalance = std::abs(change); + } + // We dampen the nodal pressure change during one iteration since our nodal pressure calculation + // is somewhat explicit. There is a relative dampening factor applied to the update value, and also + // the maximum update is limited (to 5 bar by default, can be changed with --network-max-pressure-update-in-bars). + const Scalar damped_change = std::min(damping_factor * std::abs(change), upper_update_bound); + const Scalar sign = change > 0 ? 1. : -1.; + new_pressure = pressure + sign * damped_change; } - // We dampen the nodal pressure change during one iteration since our nodal pressure calculation - // is somewhat explicit. There is a relative dampening factor applied to the update value, and also - // the maximum update is limited (to 5 bar by default, can be changed with --network-max-pressure-update-in-bars). - const Scalar damped_change = std::min(damping_factor * std::abs(change), upper_update_bound); - const Scalar sign = change > 0 ? 1. : -1.; - node_pressures_[name] = pressure + sign * damped_change; + continue; } - } else { - for (const auto& [name, pressure]: node_pressures_) { + + for (const auto& [name, pressure]: domain_pressures) { if (std::abs(pressure) > network_imbalance) { network_imbalance = std::abs(pressure); } } } + this->syncLegacyProductionState_(); for (auto& well : well_model_.genericWells()) { - - // Producers only, since we so far only support the - // "extended" network model (properties defined by - // BRANPROP and NODEPROP) which only applies to producers. + if (!well->wellEcl().predictionMode()) { + continue; + } + const auto domain = details::domainForWell(*well); + if (!domain.has_value()) { + continue; + } // The network cannot put a well without a VFP table under THP // control, so no THP limit is imposed on such a well, while its // rates still contribute to the network flows. - if (well->isProducer() && well->wellEcl().predictionMode() - && well->wellEcl().vfp_table_number() > 0) { - const auto it = node_pressures_.find(well->wellEcl().groupName()); - if (it != node_pressures_.end()) { + if (well->wellEcl().vfp_table_number() > 0) { + const auto it = this->nodePressures(*domain).find(well->wellEcl().groupName()); + if (it != this->nodePressures(*domain).end()) { + if (this->invalidNodes(*domain).count(well->wellEcl().groupName()) > 0) { + // No valid leaf pressure this iteration; keep the well's current THP limit. + continue; + } // The well belongs to a group that has a network pressure constraint; // set the dynamic THP constraint of the well accordingly. - this->imposeWellThpLimit(*well, it->second); - SingleWellState& ws = well_model_.wellState()[well->indexOfWell()]; - const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP; - // TODO: not sure why the thp is NOT updated properly elsewhere - if (thp_is_limit) { - ws.thp = well->getTHPConstraint(well_model_.summaryState()); + if (well->isProducer()) { + // For producers the leaf-node pressure is the group wellhead THP; + // apply it directly as a dynamic THP constraint. + this->imposeWellThpLimit(*well, it->second); + SingleWellState& ws = well_model_.wellState()[well->indexOfWell()]; + const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP; + // TODO: not sure why the thp is NOT updated properly elsewhere + if (thp_is_limit) { + ws.thp = well->getTHPConstraint(well_model_.summaryState()); + } + } else if (well->isInjector() && well->wellEcl().vfp_table_number() > 0) { + // For injectors, apply the network leaf-node pressure as a dynamic THP only + // if it falls within the individual well's VFPINJ table THP range. + // If P_leaf is outside the range, computeBhpAtThpLimitInj would extrapolate + // to invalid values and mark the well inoperable, causing a rate collapse. + const auto& inj_vfp = *well_model_.getVFPProperties().getInj(); + const int table_id = well->wellEcl().injectionControls( + well_model_.summaryState()).vfp_table_number; + if (inj_vfp.hasTable(table_id)) { + const auto& thp_axis = inj_vfp.getTable(table_id).getTHPAxis(); + const Scalar min_thp = static_cast(thp_axis.front()); + const Scalar max_thp = static_cast(thp_axis.back()); + const Scalar new_limit = it->second; + if (new_limit >= min_thp && new_limit <= max_thp) { + this->imposeWellThpLimit(*well, new_limit); + SingleWellState& ws = + well_model_.wellState()[well->indexOfWell()]; + const bool thp_is_limit = + ws.injection_cmode == Well::InjectorCMode::THP; + if (thp_is_limit) { + ws.thp = well->getTHPConstraint(well_model_.summaryState()); + } + } + } } } } @@ -292,12 +446,13 @@ assignNodeAndBranchValues(std::map& nodevalues, return; } - auto converged_pressures = node_pressures_; - std::tie(converged_pressures, converged_branchvalues) = this->computePressures(network, - *well_model_.getVFPProperties().getProd(), - well_model_.schedule().getUnits(), - reportStepIdx, - well_model_.comm()); + auto converged = this->computePressures(network, + *well_model_.getVFPProperties().getProd(), + well_model_.schedule().getUnits(), + reportStepIdx, + well_model_.comm()); + const auto& converged_pressures = converged.node_pressures; + converged_branchvalues = std::move(converged.branch_data); for (const auto& [node, converged_pressure] : converged_pressures) { auto it = nodevalues.find(node); assert(it != nodevalues.end() ); @@ -320,24 +475,42 @@ template void BlackoilWellModelNetworkGeneric:: initialize(const int report_step) { - // Discard pressures for nodes that are absent from the current network. - // Retained per-well limits are kept because they can outlive the network. - const auto& network = well_model_.schedule()[report_step].network(); - if (!network.active()) { - this->node_pressures_.clear(); - this->last_valid_node_pressures_.clear(); - } - else { - const auto is_stale = [&network](const auto& node_pressure) - { return !network.has_node(node_pressure.first); }; - - std::erase_if(this->node_pressures_, is_stale); - std::erase_if(this->last_valid_node_pressures_, is_stale); + const auto networks = details::networks(well_model_.schedule(), report_step); + for (const auto& [domain, network] : networks) { + // Discard pressures for nodes that are absent from the current network. + // Retained per-well limits are kept because they can outlive the network. + auto& node_pressures = this->nodePressures(domain); + auto& branch_data = this->branchData(domain); + auto& invalid_nodes = this->invalidNodes(domain); + auto& last_valid_node_pressures = this->last_valid_domain_node_pressures_[details::domainIndex(domain)]; + auto& last_valid_branch_data = this->last_valid_domain_branch_data_[details::domainIndex(domain)]; + invalid_nodes.clear(); + if (!network.get().active()) { + node_pressures.clear(); + branch_data.clear(); + last_valid_node_pressures.clear(); + last_valid_branch_data.clear(); + } + else { + const auto is_stale = [&network](const auto& node_pressure) + { return !network.get().has_node(node_pressure.first); }; + + std::erase_if(node_pressures, is_stale); + std::erase_if(branch_data, is_stale); + std::erase_if(last_valid_node_pressures, is_stale); + std::erase_if(last_valid_branch_data, is_stale); + std::erase_if(invalid_nodes, [&network](const auto& node) + { return !network.get().has_node(node); }); + } + if (domain == details::NetworkDomain::Production) { + this->syncLegacyProductionState_(); + } } + this->invalid_nodes_report_step_ = -1; // Retained THP limits can outlive network activity, so initialize every well. for (auto& well : well_model_.genericWells()) { - initializeWell(*well); + this->initializeWell(*well); } } @@ -345,22 +518,26 @@ template void BlackoilWellModelNetworkGeneric:: initializeWell(WellInterfaceGeneric& well) { - // Extended networks defined by BRANPROP and NODEPROP currently apply only - // to producers. The network cannot put a well without a VFP table under + // The network cannot put a well without a VFP table under // THP control, so no THP limit is imposed on or retained for such a well, // while its rates still contribute to the network flows. - if (!well.isProducer() || well.wellEcl().vfp_table_number() <= 0) { + if (well.wellEcl().vfp_table_number() <= 0) { return; } - const auto it = this->node_pressures_.find(well.wellEcl().groupName()); - if (it != this->node_pressures_.end()) { - // Apply and retain the network node pressure as the dynamic THP limit. - this->imposeWellThpLimit(well, it->second); - } else { - // Reapply a retained limit after network detachment or well reconstruction. - const auto& ws = well_model_.wellState().well(well.indexOfWell()); - if (ws.network_thp_limit.has_value()) { - well.setDynamicThpLimit(*ws.network_thp_limit); + + const auto domain = details::domainForWell(well); + if (domain.has_value() && !this->nodePressures(*domain).empty()) { + const auto it = this->nodePressures(*domain).find(well.wellEcl().groupName()); + if (it != this->nodePressures(*domain).end() && well.isProducer()) { + // Carry the converged production-network pressure into the next + // report step as the producer's starting THP constraint. + this->imposeWellThpLimit(well, it->second); + } else if (it == this->nodePressures(*domain).end()) { + // Reapply a retained limit after network detachment or well reconstruction. + const auto& ws = well_model_.wellState().well(well.indexOfWell()); + if (ws.network_thp_limit.has_value()) { + well.setDynamicThpLimit(*ws.network_thp_limit); + } } } } @@ -374,7 +551,7 @@ imposeWellThpLimit(WellInterfaceGeneric& well, const Scalar } template -std::pair, std::map> +typename BlackoilWellModelNetworkGeneric::NetworkPressures BlackoilWellModelNetworkGeneric:: computePressures(const Network::ExtNetwork& network, const VFPProdProperties& vfp_prod_props, @@ -392,7 +569,31 @@ computePressures(const Network::ExtNetwork& network, network_pressure_computation( well_model_, network, vfp_prod_props, unit_system, reportStepIdx, comm); - return network_pressure_computation.run(); + auto [node_pressures, branch_data] = network_pressure_computation.run(); + return {std::move(node_pressures), std::move(branch_data), network_pressure_computation.invalidNodes()}; +} + +template +typename BlackoilWellModelNetworkGeneric::NetworkPressures +BlackoilWellModelNetworkGeneric:: +computePressures(const Network::ExtNetwork& network, + const VFPInjProperties& vfp_inj_props, + const UnitSystem& unit_system, + const int reportStepIdx, + const Parallel::Communication& comm) const +{ + OPM_TIMEFUNCTION(); + if (!network.active()) { + return {}; + } + + NetworkPressureComputation, + VFPInjProperties> + network_pressure_computation( + well_model_, network, vfp_inj_props, unit_system, reportStepIdx, comm); + + auto [node_pressures, branch_data] = network_pressure_computation.run(); + return {std::move(node_pressures), std::move(branch_data), network_pressure_computation.invalidNodes()}; } template @@ -404,7 +605,11 @@ operator==(const BlackoilWellModelNetworkGeneric& rhs) const && this->node_pressures_ == rhs.node_pressures_ && this->last_valid_node_pressures_ == rhs.last_valid_node_pressures_ && this->branch_data_ == rhs.branch_data_ - && this->last_valid_branch_data_ == rhs.last_valid_branch_data_; + && this->last_valid_branch_data_ == rhs.last_valid_branch_data_ + && this->domain_node_pressures_ == rhs.domain_node_pressures_ + && this->last_valid_domain_node_pressures_ == rhs.last_valid_domain_node_pressures_ + && this->domain_branch_data_ == rhs.domain_branch_data_ + && this->last_valid_domain_branch_data_ == rhs.last_valid_domain_branch_data_; } template class BlackoilWellModelNetworkGeneric; diff --git a/opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp index 39a9a3ef149..f044fc2bd4a 100644 --- a/opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp @@ -24,14 +24,17 @@ #define OPM_BLACKOILWELLMODEL_NETWORK_GENERIC_HEADER_INCLUDED #include +#include #include #include #include +#include #include #include +#include #include namespace Opm { @@ -39,11 +42,65 @@ namespace Opm { class UnitSystem; template class BlackoilWellModelGeneric; template class WellInterfaceGeneric; + template class VFPInjProperties; template class VFPProdProperties; } namespace Opm { +namespace details { + + enum class NetworkDomain : std::size_t { + Production = 0, + InjectionGas, + InjectionWater, + Count + }; + + constexpr std::size_t domainIndex(const NetworkDomain domain) + { + return static_cast(domain); + } + + struct NetworkDescriptor { + NetworkDomain domain; + std::reference_wrapper network; + }; + + /// The network domain corresponding to a well's type + /// (producer, or injector and water/gas injection phase), or nullopt. + template + std::optional domainForWell(const Well& well) + { + if (well.isProducer()) { + return NetworkDomain::Production; + } + if (well.isInjector()) { + if (well.wellEcl().injectorType() == InjectorType::GAS) { + return NetworkDomain::InjectionGas; + } + if (well.wellEcl().injectorType() == InjectorType::WATER) { + return NetworkDomain::InjectionWater; + } + } + return std::nullopt; + } + + /// Helper to check if any network (production, gas injection, water injection) is active at a given time step. + bool anyNetworkActive(const Schedule& schedule, const int timeStepIdx); + + /// Helper to get all active networks (production, gas injection, water injection) at a given time step. + std::vector + activeNetworks(const Schedule& schedule, const int timeStepIdx); + + /// Helper to get all networks (production, gas injection, water injection) at a given time step, + /// whether active or not. + std::vector + networks(const Schedule& schedule, const int timeStepIdx); + +} // namespace details + + /// Class for handling the blackoil well network model. template class BlackoilWellModelNetworkGeneric @@ -99,12 +156,16 @@ class BlackoilWellModelNetworkGeneric { this->last_valid_node_pressures_ = this->node_pressures_; this->last_valid_branch_data_ = this->branch_data_; + this->last_valid_domain_node_pressures_ = this->domain_node_pressures_; + this->last_valid_domain_branch_data_ = this->domain_branch_data_; } void resetState() { this->node_pressures_ = this->last_valid_node_pressures_; this->branch_data_ = this->last_valid_branch_data_; + this->domain_node_pressures_ = this->last_valid_domain_node_pressures_; + this->domain_branch_data_ = this->last_valid_domain_branch_data_; } template @@ -114,23 +175,92 @@ class BlackoilWellModelNetworkGeneric serializer(last_valid_node_pressures_); serializer(branch_data_); serializer(last_valid_branch_data_); + serializer(domain_node_pressures_); + serializer(last_valid_domain_node_pressures_); + serializer(domain_branch_data_); + serializer(last_valid_domain_branch_data_); } bool operator==(const BlackoilWellModelNetworkGeneric& rhs) const; protected: - std::pair, std::map> + //! \brief Apply a network THP limit and retain it across well + //! reconstruction or network detachment. + void imposeWellThpLimit(WellInterfaceGeneric& well, + const Scalar limit); + + /// Result of one network pressure evaluation for one network (domain). + struct NetworkPressures + { + std::map node_pressures; + std::map branch_data; + // Nodes (and their descendants) whose VFP lookup has no solution; their + // node_pressures entries are placeholders and must not be used. + std::set invalid_nodes; + }; + + NetworkPressures computePressures(const Network::ExtNetwork& network, const VFPProdProperties& vfp_prod_props, const UnitSystem& unit_system, const int reportStepIdx, const Parallel::Communication& comm) const; - //! \brief Apply a network THP limit and retain it across well - //! reconstruction or network detachment. - void imposeWellThpLimit(WellInterfaceGeneric& well, - const Scalar limit); + NetworkPressures + computePressures(const Network::ExtNetwork& network, + const VFPInjProperties& vfp_inj_props, + const UnitSystem& unit_system, + const int reportStepIdx, + const Parallel::Communication& comm) const; + void updateActiveStateImpl(const Network::ExtNetwork& network); + + static constexpr details::NetworkDomain productionNetworkDomain() + { + return details::NetworkDomain::Production; + } + + const std::map& nodePressures(const details::NetworkDomain domain) const + { + return domain_node_pressures_[details::domainIndex(domain)]; + } + + std::map& nodePressures(const details::NetworkDomain domain) + { + return domain_node_pressures_[details::domainIndex(domain)]; + } + + const std::map& branchData(const details::NetworkDomain domain) const + { + return domain_branch_data_[details::domainIndex(domain)]; + } + + std::map& branchData(const details::NetworkDomain domain) + { + return domain_branch_data_[details::domainIndex(domain)]; + } + + const std::set& invalidNodes(const details::NetworkDomain domain) const + { + return domain_invalid_nodes_[details::domainIndex(domain)]; + } + + std::set& invalidNodes(const details::NetworkDomain domain) + { + return domain_invalid_nodes_[details::domainIndex(domain)]; + } + + void syncLegacyProductionState_() + { + this->node_pressures_ = this->nodePressures(productionNetworkDomain()); + this->branch_data_ = this->branchData(productionNetworkDomain()); + } + + void syncProductionDomainState_() + { + this->nodePressures(productionNetworkDomain()) = this->node_pressures_; + this->branchData(productionNetworkDomain()) = this->branch_data_; + } bool active_{false}; BlackoilWellModelGeneric& well_model_; @@ -139,10 +269,19 @@ class BlackoilWellModelNetworkGeneric std::map node_pressures_; // Network branch pressure drops and flow rates for output (outlet branch for production network, inlet branch for injection network) std::map branch_data_; + // Domain-scoped pressure state to avoid collisions between production and injection networks. + std::array, details::domainIndex(details::NetworkDomain::Count)> domain_node_pressures_; + std::array, details::domainIndex(details::NetworkDomain::Count)> domain_branch_data_; + // Nodes without a valid VFP solution in the last evaluation (per domain); not serialized, + // recomputed on every updatePressures(). + std::array, details::domainIndex(details::NetworkDomain::Count)> domain_invalid_nodes_; + int invalid_nodes_report_step_{-1}; // Valid network pressures for output and initialization for safe restart after failed iterations std::map last_valid_node_pressures_; // Valid network branch pressure drops and flow rates for output (outlet branch for production network, inlet branch for injection network) for safe restart after failed iterations std::map last_valid_branch_data_; + std::array, details::domainIndex(details::NetworkDomain::Count)> last_valid_domain_node_pressures_; + std::array, details::domainIndex(details::NetworkDomain::Count)> last_valid_domain_branch_data_; }; } // namespace Opm diff --git a/opm/simulators/wells/BlackoilWellModelNetworkPressureComputation.hpp b/opm/simulators/wells/BlackoilWellModelNetworkPressureComputation.hpp index a63dbe175f5..5d59d1757c7 100644 --- a/opm/simulators/wells/BlackoilWellModelNetworkPressureComputation.hpp +++ b/opm/simulators/wells/BlackoilWellModelNetworkPressureComputation.hpp @@ -25,15 +25,20 @@ #include #include #include +#include #include #include +#include #include #include +#include + #include #include +#include #include #include #include @@ -43,6 +48,48 @@ namespace Opm { +/// Result of a single network branch VFP lookup. +template +struct NetworkBranchPressure +{ + Scalar pressure{0.0}; + // False when the table has no solution at this point (zero-filled cells give + // bhp <= 1 atm); the pressure must then not be used as a node pressure. + bool valid{true}; + // True when the flow rate or upstream pressure had to be clamped to the table axes. + bool clamped{false}; +}; + +namespace detail { + /// Clamp the VFP lookup point to the table axes; the tables must not be extrapolated + /// for network branches (a zero-filled tail extrapolates to negative pressures). + /// Rates are scaled uniformly so that WFR/GFR fractions are preserved. + template + bool clampToTableAxes(const Table& table, std::vector& rates, Scalar& up_press) + { + bool clamped = false; + const auto& thp_axis = table.getTHPAxis(); + const Scalar thp_lo = thp_axis.front(); + const Scalar thp_hi = thp_axis.back(); + if (up_press < thp_lo || up_press > thp_hi) { + up_press = std::clamp(up_press, thp_lo, thp_hi); + clamped = true; + } + const auto& flo_axis = table.getFloAxis(); + const Scalar flo = std::abs(getFlo(table, + rates[IndexTraits::waterPhaseIdx], + rates[IndexTraits::oilPhaseIdx], + rates[IndexTraits::gasPhaseIdx])); + const Scalar flo_hi = flo_axis.back(); + if (flo > flo_hi && flo > 0.0) { + const Scalar s = flo_hi / flo; + std::ranges::transform(rates, rates.begin(), [s](const auto r) { return s * r; }); + clamped = true; + } + return clamped; + } +} // namespace detail + /// @brief Helper class to insulate the NetworkPressureComputation class from /// the differences between production and injection VFP tables. template @@ -58,37 +105,53 @@ struct NetworkVfpPressureCalculator + static bool hasLeafNodeRate(const GroupState& group_state, + const std::string& node) + { + return group_state.has_network_leaf_node_production_rates(node); + } + template static const std::vector - leafNodeRate(const GroupState& group_state, const std::string& node) + leafNodeRate(const GroupState& group_state, + const std::string& node) { return group_state.network_leaf_node_production_rates(node); } template - static Scalar compute(const VFPProdProperties& vfp_props, - const int table_id, - const std::vector& rates, - const Scalar up_press, - const Branch& upbranch, - const UnitSystem& unit_system) + static NetworkBranchPressure compute(const VFPProdProperties& vfp_props, + const int table_id, + std::vector rates, + Scalar up_press, + const Branch& upbranch, + const UnitSystem& unit_system) { // 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. - const auto alq_type = vfp_props.getTable(table_id).getALQType(); + const auto& table = vfp_props.getTable(table_id); + const auto alq_type = table.getALQType(); const auto dimension = VFPProdTable::ALQDimension(alq_type, unit_system); const Scalar alq = upbranch.alq_value(dimension).value_or(0.0); - return vfp_props.bhp(table_id, - 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 + NetworkBranchPressure result; + // Preserve the established production-network behaviour. Production VFP + // tables have historically been extrapolated outside their axes, and + // existing production cases rely on that when the result remains valid. + result.clamped = false; + result.pressure = vfp_props.bhp(table_id, + 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 + result.valid = result.pressure > unit::atm; + return result; } }; @@ -100,26 +163,38 @@ struct NetworkVfpPressureCalculator + static bool hasLeafNodeRate(const GroupState& group_state, + const std::string& node) + { + return group_state.has_network_leaf_node_injection_rates(node); + } + template static const std::vector - leafNodeRate(const GroupState& group_state, const std::string& node) + leafNodeRate(const GroupState& group_state, + const std::string& node) { return group_state.network_leaf_node_injection_rates(node); } template - static Scalar compute(const VFPInjProperties& vfp_props, - const int table_id, - const std::vector& rates, - const Scalar up_press, - const Branch&, - const UnitSystem&) + static NetworkBranchPressure compute(const VFPInjProperties& vfp_props, + const int table_id, + std::vector rates, + Scalar up_press, + const Branch&, + const UnitSystem&) { - return vfp_props.bhp(table_id, - rates[IndexTraits::waterPhaseIdx], - rates[IndexTraits::oilPhaseIdx], - rates[IndexTraits::gasPhaseIdx], - up_press); + NetworkBranchPressure result; + result.clamped = detail::clampToTableAxes(vfp_props.getTable(table_id), rates, up_press); + result.pressure = vfp_props.bhp(table_id, + rates[IndexTraits::waterPhaseIdx], + rates[IndexTraits::oilPhaseIdx], + rates[IndexTraits::gasPhaseIdx], + up_press); + result.valid = result.pressure > unit::atm; + return result; } }; @@ -171,11 +246,34 @@ class NetworkPressureComputation // Going the other way (from roots to leafs), calculate the pressure // at each node using VFP tables and rates. computeNodePressures(root_to_child_nodes, node_inflows); + +#ifdef OPM_NETWORK_PRESSURE_TRACE + // Off unless the macro is defined: this builds a string per node per + // sub-iteration per domain, whether or not the log keeps it. + OpmLog::debug("Network pressure computation completed for root " + root.get().name() + ". Node pressures:"); + for (const auto& [node, pressure] : node_pressures_) { + OpmLog::debug("Network node " + node + " pressure: " + std::to_string(pressure/1e5) + " bar"); + } + OpmLog::debug("Node inflows:"); + for (const auto& [node, inflows] : node_inflows) { + OpmLog::debug("Network node " + node + " inflows: " + + std::to_string(inflows[0]*86400) + ", " + std::to_string(inflows[1]*86400) + ", " + std::to_string(inflows[2]*86400)); + } +#endif + } return {node_pressures_, branch_data_}; } + /// Nodes whose pressure could not be computed from the VFP tables (and their + /// descendants). Their entries in the pressure map are placeholders (the upstream + /// pressure) and must not be used as node pressures. + const std::set& invalidNodes() const + { + return invalid_nodes_; + } + private: std::pair, std::set> collectTreeNodes(const std::string& root) const @@ -208,14 +306,18 @@ class NetworkPressureComputation 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)) { + // Guard against empty leaf nodes (may not be present in GRUPTREE). + // Use the domain-correct check so injection networks query the injection + // rate map rather than the production rate map (which is always empty for + // pure injection groups, causing zero-rate pressure calculations). + using Calc = NetworkVfpPressureCalculator; + if (!Calc::hasLeafNodeRate(well_model_.groupStateHelper().groupState(), node)) { node_inflows[node] = zero_rates; continue; } - using Calc = NetworkVfpPressureCalculator; - node_inflows[node] = Calc::leafNodeRate(well_model_.groupStateHelper().groupState(), node); + node_inflows[node] = Calc::leafNodeRate(well_model_.groupStateHelper().groupState(), + node); if (network_.node(node).add_gas_lift_gas()) { addGasLiftGas(node, node_inflows[node]); } @@ -248,7 +350,7 @@ class NetworkPressureComputation } // Sum ALQ across all processes to get total ALQ for the node. // Note that communication is required here since each - // process has different wells, and the loop above therefore + // process has different wells, and the loop above therefore // only considers local wells. // However, all processes have all groups and their rates available, // so we do not need to communicate those. @@ -319,7 +421,12 @@ class NetworkPressureComputation continue; } - const Scalar up_press = node_pressures_[(*upbranch).uptree_node()]; + const std::string& up_node = (*upbranch).uptree_node(); + const Scalar up_press = node_pressures_[up_node]; + // Descendants of a node without a valid pressure have none either. + if (invalid_nodes_.count(up_node) > 0) { + invalid_nodes_.insert(node); + } const auto vfp_table = (*upbranch).vfp_table(); if (!vfp_table) { // Table number specified as 9999 in the deck, no pressure loss. @@ -342,7 +449,23 @@ class NetworkPressureComputation auto rates = node_inflows.at(node); assert(rates.size() == 3); Calc::prepareRates(rates); - auto node_pressure = Calc::compute(vfp_props_, *vfp_table, rates, up_press, *upbranch, unit_system_); + const auto branch = Calc::compute(vfp_props_, *vfp_table, rates, up_press, *upbranch, unit_system_); + // An invalid lookup (zero-filled table cells) gets the upstream pressure as a + // placeholder so downstream lookups stay in range; callers must consult invalidNodes(). + const Scalar node_pressure = branch.valid ? branch.pressure : up_press; + if (!branch.valid) { + invalid_nodes_.insert(node); + } + if (!branch.valid || branch.clamped) { + OpmLog::debug(fmt::format("Network branch {} -> {}: VFP table {} {} at rates ({:.4g}, {:.4g}, {:.4g}) sm3/d, " + "upstream pressure {:.2f} bar", + up_node, node, *vfp_table, + branch.valid ? "lookup clamped to the table axes" : "has no solution", + rates[IndexTraits::waterPhaseIdx] * unit::day, + rates[IndexTraits::oilPhaseIdx] * unit::day, + rates[IndexTraits::gasPhaseIdx] * unit::day, + up_press / unit::barsa)); + } node_pressures_[node] = node_pressure; // Prefer inserting after computing the pressure, hence negating rates branch_data_.try_emplace(node, @@ -361,6 +484,7 @@ class NetworkPressureComputation const Communication& comm_; std::map node_pressures_; std::map branch_data_; + std::set invalid_nodes_; }; } // namespace Opm diff --git a/opm/simulators/wells/BlackoilWellModelNetwork_impl.hpp b/opm/simulators/wells/BlackoilWellModelNetwork_impl.hpp index 90aa5845696..6964f171274 100644 --- a/opm/simulators/wells/BlackoilWellModelNetwork_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModelNetwork_impl.hpp @@ -41,6 +41,8 @@ #include +#include + namespace Opm { template @@ -93,8 +95,7 @@ update(const bool mandatory_network_balance, { OPM_TIMEFUNCTION(); const int episodeIdx = well_model_.simulator().episodeIndex(); - const auto& network = well_model_.schedule()[episodeIdx].network(); - if (!well_model_.wellsActive() && !network.active()) { + if (!well_model_.wellsActive() && !details::anyNetworkActive(well_model_.schedule(), episodeIdx)) { return {/*more_network_update=*/false, /*network_imbalance=*/0.0}; } @@ -140,16 +141,41 @@ update(const bool mandatory_network_balance, } for (const auto& well : well_model_) { - if (well->isInjector() || !well->wellEcl().predictionMode()) { + if (!well->wellEcl().predictionMode()) { continue; } - const auto it = this->node_pressures_.find(well->wellEcl().groupName()); - if (it != this->node_pressures_.end()) { + const auto domain = details::domainForWell(*well); + + if (!domain.has_value()) { + continue; + } + + const auto it = this->nodePressures(*domain).find(well->wellEcl().groupName()); + if (it != this->nodePressures(*domain).end()) { well->prepareWellBeforeAssembling(well_model_.simulator(), dt, well_model_.groupStateHelper(), well_model_.wellState()); + // Option B: after re-solving at the current network THP, update + // ws.well_potentials for injection wells. The rate_less_than_potential + // check in WellConstraints::activeInjectionConstraint compares current + // injection rates against ws.well_potentials to decide whether switching + // to THP mode would increase or decrease injection. The potentials are + // normally computed once per timestep at the static WCONINJE THP and are + // stale during network iterations. Refreshing them here (from the rate + // the well just solved to under the current network THP) makes the check + // accurate for subsequent outer iterations. + if (well->isInjector()) { + auto& ws = well_model_.wellState().well(well->indexOfWell()); + if (ws.injection_cmode == Well::InjectorCMode::THP) { + const int np = well_model_.numPhases(); + for (int p = 0; p < np; ++p) { + ws.well_potentials[p] = + std::max(Scalar{0.0}, ws.surface_rates[p]); + } + } + } } } well_model_.updateAndCommunicateGroupData(episodeIdx, /*update_wellgrouptarget*/ true); @@ -166,6 +192,10 @@ computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger) { OPM_TIMEFUNCTION(); const int reportStepIdx = well_model_.simulator().episodeIndex(); + // This function is only relevant for auto-choke groups, and + // therefore as of now only relevant for the production network. + // \TODO: If we later also want to support auto-choke groups in the + // injection network, we should change this function also. const auto& network = well_model_.schedule()[reportStepIdx].network(); const auto& balance = well_model_.schedule()[reportStepIdx].network_balance(); const Scalar thp_tolerance = balance.thp_tolerance(); diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 5fa224b9967..942eb04efa6 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -603,6 +603,7 @@ namespace Opm { } this->network_.initializeWell(*well); + try { using GLiftEclWells = typename GasLiftGroupInfo::GLiftEclWells; GLiftEclWells ecl_well_map; @@ -1175,8 +1176,7 @@ namespace Opm { this->updateNetworkActiveState_(); } const int episodeIdx = simulator_.episodeIndex(); - const auto& network = this->schedule()[episodeIdx].network(); - if (!this->wellsActive() && !network.active()) { + if (!this->wellsActive() && !details::anyNetworkActive(this->schedule(), episodeIdx)) { return; } } @@ -1243,8 +1243,15 @@ namespace Opm { while (do_network_update) { if (!this->isRescoupSlaveCoupledNetworkIteration_() && network_update_iteration >= max_iteration ) { - // only output to terminal if we at the last newton iterations where we try to balance the network. const int episodeIdx = simulator_.episodeIndex(); + const auto& balance = this->schedule()[episodeIdx].network_balance(); + // If the imbalance is already within tolerance, the network is converged; the + // outer loop continued only due to ALQ or control changes. Don't report this + // as an unconverged result -- just break quietly. + if (network_imbalance <= balance.pressure_tolerance()) { + break; + } + // only output to terminal if we at the last newton iterations where we try to balance the network. if (this->network_.willBalanceOnNextIteration(episodeIdx)) { if (this->terminal_output_) { const std::string msg = fmt::format("Maximum of {:d} network iterations has been used and we stop the update, \n" diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index c0745df5fc6..5fb7418b004 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -102,6 +102,12 @@ void GroupState::update_production_rates(const std::string& gname, this->m_production_rates[gname] = rates; } +template +bool GroupState::has_network_leaf_node_injection_rates(const std::string& gname) const +{ + return this->m_network_leaf_node_injection_rates.count(gname) > 0; +} + template void GroupState::update_network_leaf_node_injection_rates(const std::string& gname, const std::vector& rates) @@ -165,6 +171,12 @@ GroupState::network_leaf_node_injection_rates(const std::string& gname) return group_iter->second; } +template +bool GroupState::has_network_leaf_node_production_rates(const std::string& gname) const +{ + return this->m_network_leaf_node_production_rates.count(gname) > 0; +} + template const std::vector& GroupState::network_leaf_node_production_rates(const std::string& gname) const diff --git a/opm/simulators/wells/GroupState.hpp b/opm/simulators/wells/GroupState.hpp index e7267af84de..c0a4110a9a8 100644 --- a/opm/simulators/wells/GroupState.hpp +++ b/opm/simulators/wells/GroupState.hpp @@ -55,7 +55,9 @@ class GroupState { void update_network_leaf_node_production_rates(const std::string& gname, const std::vector& rates); const std::vector& production_rates(const std::string& gname) const; + bool has_network_leaf_node_injection_rates(const std::string& gname) const; const std::vector& network_leaf_node_injection_rates(const std::string& gname) const; + bool has_network_leaf_node_production_rates(const std::string& gname) const; const std::vector& network_leaf_node_production_rates(const std::string& gname) const; void update_well_group_thp(const std::string& gname, const double& thp); diff --git a/opm/simulators/wells/GroupStateHelper.cpp b/opm/simulators/wells/GroupStateHelper.cpp index 67c086fa5a2..28570735259 100644 --- a/opm/simulators/wells/GroupStateHelper.cpp +++ b/opm/simulators/wells/GroupStateHelper.cpp @@ -1280,9 +1280,13 @@ GroupStateHelper::updateNetworkLeafNodeRates() } }; do_update(this->schedule_[this->report_step_].network(), /*is_injector=*/false); - // TODO: do the below to support injection networks when available. - // do_update(this->schedule_[this->report_step_].gas_injection_network(), /*is_injector=*/true); - // do_update(this->schedule_[this->report_step_].water_injection_network(), /*is_injector=*/true); + for (const Phase phase : {Phase::GAS, Phase::WATER}) { + if (const auto injNetwork = this->schedule_[this->report_step_].injectionNetwork.get_ptr(phase); + injNetwork != nullptr) + { + do_update(*injNetwork, /* is_injector = */ true); + } + } } template diff --git a/opm/simulators/wells/WellConstraints.cpp b/opm/simulators/wells/WellConstraints.cpp index c358c1d1179..57d01a025ed 100644 --- a/opm/simulators/wells/WellConstraints.cpp +++ b/opm/simulators/wells/WellConstraints.cpp @@ -148,20 +148,30 @@ activeInjectionConstraint(const SingleWellState& ws, return Well::InjectorCMode::RESV; } - // Note: we are not working on injecting network yet, so it is possible we need to change the following line - // to be as follows to incorporate the injecting network nodal pressure - // if (well_.wellHasTHPConstraints(summaryState) && currentControl != Well::InjectorCMode::THP) - if (controls.hasControl(Well::InjectorCMode::THP) && currentControl != Well::InjectorCMode::THP) + // Use wellHasTHPConstraints so that injection wells with a dynamic THP from + // the injection network (dynamic_thp_limit_ set) also enter this check. + // Wells with neither an explicit WCONINJE THP nor a network-derived THP are + // unaffected because wellHasTHPConstraints returns false for them. + if (well_.wellHasTHPConstraints(summaryState) && currentControl != Well::InjectorCMode::THP) { const auto& thp = well_.getTHPConstraint(summaryState); Scalar current_thp = ws.thp; if (thp < current_thp) { + // When the THP comes from the injection network (dynamic_thp_limit_ is set), + // well potentials were computed before network iterations at a different (static) + // THP and are stale. The rate_less_than_potential check would always suppress + // switching in that case. Bypass the check for dynamic THP — this mirrors the + // default production-well behaviour (no WVFPEXP) where switching is unconditional. bool rate_less_than_potential = true; - for (int p = 0; p < well_.numPhases(); ++p) { - // Currently we use the well potentials here computed before the iterations. - // We may need to recompute the well potentials to get a more - // accurate check here. - rate_less_than_potential = rate_less_than_potential && (ws.surface_rates[p]) <= ws.well_potentials[p]; + if (!well_.getDynamicThpLimit().has_value()) { + for (int p = 0; p < well_.numPhases(); ++p) { + // Currently we use the well potentials here computed before the iterations. + // We may need to recompute the well potentials to get a more + // accurate check here. + rate_less_than_potential = rate_less_than_potential && (ws.surface_rates[p]) <= ws.well_potentials[p]; + } + } else { + rate_less_than_potential = false; } if (!rate_less_than_potential) { thp_limit_violated_but_not_switched = false; diff --git a/tests/test_networkpressure.cpp b/tests/test_networkpressure.cpp index a2ee76ced92..79b130512bf 100644 --- a/tests/test_networkpressure.cpp +++ b/tests/test_networkpressure.cpp @@ -44,6 +44,7 @@ #include #include #include +#include #include #include #include @@ -228,20 +229,28 @@ struct MockWellModel struct MockGroupState { + // Leaf rates in Sm3/day, phase order water, oil, gas. Tests may override + // these before running the computation. + static inline std::vector injection_rates_sm3_day {500.0, 0.0, 5000.0}; + static inline std::vector production_rates_sm3_day {500.0, 500.0, 5000.0}; + + static std::vector toSI(const std::vector& r) + { + std::vector out(r.size()); + std::ranges::transform(r, out.begin(), [](double v) { return convert::from(v, cubic(meter) / day); }); + return out; + } + bool has_production_rates(const std::string) const { return true; } + bool has_network_leaf_node_injection_rates(const std::string) const { return true; } + bool has_network_leaf_node_production_rates(const std::string) const { return true; } std::vector network_leaf_node_injection_rates(const std::string) const { - // Phase order water, oil, gas. - return {convert::from(500.0, cubic(meter) / day), - 0.0, - convert::from(5000.0, cubic(meter) / day)}; + return toSI(injection_rates_sm3_day); } std::vector network_leaf_node_production_rates(const std::string) const { - // Phase order water, oil, gas. - return {convert::from(500.0, cubic(meter) / day), - convert::from(500.0, cubic(meter) / day), - convert::from(5000.0, cubic(meter) / day)}; + return toSI(production_rates_sm3_day); } Scalar well_group_thp(const std::string&) const { return convert::from(100.0, bars); } }; @@ -290,7 +299,7 @@ double terminalPressure(NetworkScenario scenario) struct NetworkSetup { - NetworkSetup(NetworkScenario scenario) + NetworkSetup(NetworkScenario scenario, std::optional terminal_pressure_override = std::nullopt) : deck{Parser{}.parseString(inputString(scenario))} { // Set up VFP property objects. @@ -306,8 +315,11 @@ struct NetworkSetup network.add_branch(Network::Branch{"M5S", "PLAT-A", 3, 0.0}); network.add_branch(Network::Branch{"G1", "M5S", 9999, 0.0}); Network::Node node{"PLAT-A"}; - node.terminal_pressure(terminalPressure(scenario)); + node.terminal_pressure(terminal_pressure_override.value_or(terminalPressure(scenario))); network.update_node(node); + // Restore the default leaf rates so tests do not leak state into each other. + MockWellModel::MockGroupStateHelper::MockGroupState::injection_rates_sm3_day = {500.0, 0.0, 5000.0}; + MockWellModel::MockGroupStateHelper::MockGroupState::production_rates_sm3_day = {500.0, 500.0, 5000.0}; } Deck deck; @@ -332,7 +344,7 @@ BOOST_AUTO_TEST_CASE(gas_injection_pressure_computation) BOOST_CHECK_CLOSE(s.vfp_inj_props.bhp(3, 0.0, 0.0, gasrate, thp), expected_bhp, 1e-7); using Comm = Dune::Communication; - // NetworkPressureComputation stores const references to comm and unit system, hence + // NetworkPressureComputation stores const references to comm and unit system, hence // we need to make sure that their lifetime is longer than the constructor lasts auto comm = Comm{}; auto unit_system = UnitSystem {}; @@ -357,7 +369,7 @@ BOOST_AUTO_TEST_CASE(water_injection_pressure_computation) // Test using mock setup. using Comm = Dune::Communication; - // NetworkPressureComputation stores const references to comm and unit system, hence + // NetworkPressureComputation stores const references to comm and unit system, hence // we need to make sure that their lifetime is longer than the constructor lasts auto comm = Comm{}; auto unit_system = UnitSystem {}; @@ -393,4 +405,86 @@ BOOST_AUTO_TEST_CASE(production_pressure_computation) BOOST_CHECK_CLOSE(pressures.at("G1"), expected_pressure, 1e-7); } +BOOST_AUTO_TEST_CASE(production_rate_beyond_flow_axis_is_extrapolated) +{ + auto s = NetworkSetup{NetworkScenario::Production}; + // The last flow-axis value is 2000 Sm3/d. At 2500 Sm3/d and 20 bar + // THP, legacy linear extrapolation gives 45 bar; clamping gives 40 bar. + MockWellModel::MockGroupStateHelper::MockGroupState::production_rates_sm3_day = + {0.0, 2500.0, 0.0}; + + using Comm = Dune::Communication; + auto comm = Comm{}; + auto unit_system = UnitSystem{}; + NetworkPressureComputation, Comm> comp( + s.well_model, s.network, s.vfp_prod_props, unit_system, 0, comm); + const auto [pressures, branch_data] = comp.run(); + BOOST_REQUIRE(pressures.find("G1") != pressures.end()); + BOOST_CHECK_CLOSE(pressures.at("G1"), convert::from(45.0, bars), 1e-7); + BOOST_CHECK(comp.invalidNodes().empty()); +} + +// The tables below use zero-filled cells for (rate, THP) combinations the flow line +// cannot deliver, and their axes do not cover every state the wells may be in during +// network iterations. A network branch lookup must never extrapolate into that region +// (it gives negative pressures) nor accept a zero-filled cell as a node pressure. + +BOOST_AUTO_TEST_CASE(gas_injection_rate_beyond_flow_axis) +{ + auto s = NetworkSetup{NetworkScenario::GasInjection}; + // 2.5e6 Sm3/d is beyond the last flow-axis point (2.0e6); a linear extrapolation of the + // THP=350 row (..., 86.011, 0.000) gives a pressure of about -213 bar. + MockWellModel::MockGroupStateHelper::MockGroupState::injection_rates_sm3_day = {0.0, 0.0, 2.5e6}; + + using Comm = Dune::Communication; + auto comm = Comm{}; + auto unit_system = UnitSystem {}; + NetworkPressureComputation, Comm> comp( + s.well_model, s.network, s.vfp_inj_props, unit_system, 0, comm); + const auto [pressures, branch_data] = comp.run(); + BOOST_REQUIRE(pressures.find("G1") != pressures.end()); + // Clamped to the axis end the table gives 0.0 -> no solution: the node is flagged and the + // placeholder pressure is the upstream (terminal) pressure, never a negative value. + BOOST_CHECK(pressures.at("M5S") >= unit::atm); + BOOST_CHECK(pressures.at("G1") >= unit::atm); + BOOST_CHECK(comp.invalidNodes().count("M5S") == 1); + BOOST_CHECK(comp.invalidNodes().count("G1") == 1); +} + +BOOST_AUTO_TEST_CASE(gas_injection_zero_cell_region) +{ + // At THP=100 bar the table is zero for rates >= 589394 Sm3/d. + auto s = NetworkSetup{NetworkScenario::GasInjection, convert::from(100.0, bars)}; + MockWellModel::MockGroupStateHelper::MockGroupState::injection_rates_sm3_day = {0.0, 0.0, 6.0e5}; + + using Comm = Dune::Communication; + auto comm = Comm{}; + auto unit_system = UnitSystem {}; + NetworkPressureComputation, Comm> comp( + s.well_model, s.network, s.vfp_inj_props, unit_system, 0, comm); + const auto [pressures, branch_data] = comp.run(); + BOOST_REQUIRE(pressures.find("G1") != pressures.end()); + BOOST_CHECK(pressures.at("G1") >= unit::atm); + BOOST_CHECK(comp.invalidNodes().count("M5S") == 1); + BOOST_CHECK(comp.invalidNodes().count("G1") == 1); +} + +BOOST_AUTO_TEST_CASE(gas_injection_thp_below_axis) +{ + // Terminal pressure 20 bar is below the first THP-axis point (50 bar). Extrapolating the + // first interval gives 68.834 - 0.6*(135.406 - 68.834) = 28.9 bar; clamping to the axis + // gives the THP=50 row value 68.834 bar. + auto s = NetworkSetup{NetworkScenario::GasInjection, convert::from(20.0, bars)}; + + using Comm = Dune::Communication; + auto comm = Comm{}; + auto unit_system = UnitSystem {}; + NetworkPressureComputation, Comm> comp( + s.well_model, s.network, s.vfp_inj_props, unit_system, 0, comm); + const auto [pressures, branch_data] = comp.run(); + BOOST_REQUIRE(pressures.find("G1") != pressures.end()); + BOOST_CHECK_CLOSE(pressures.at("G1"), convert::from(68.834, bars), 1e-7); + BOOST_CHECK(comp.invalidNodes().empty()); +} + BOOST_AUTO_TEST_SUITE_END() // NetworkPressureComputationTests