diff --git a/opm/input/eclipse/Schedule/Network/Balance.hpp b/opm/input/eclipse/Schedule/Network/Balance.hpp index fb3120855..05b69616d 100644 --- a/opm/input/eclipse/Schedule/Network/Balance.hpp +++ b/opm/input/eclipse/Schedule/Network/Balance.hpp @@ -17,33 +17,34 @@ along with OPM. If not, see . */ - #ifndef NETWORK_BALANCE_HPP #define NETWORK_BALANCE_HPP -#include -#include +#include +#include namespace Opm { -class DeckKeyword; -struct Tuning; -class UnitSystem; + class DeckKeyword; + class UnitSystem; +} // namespace Opm -namespace Network { +namespace Opm { namespace Network { -class Balance { +class Balance +{ public: + enum class CalcMode { + None = 0, + TimeInterval = 1, + TimeStepStart = 2, + NUPCOL = 3, + }; -enum class CalcMode { - None = 0, - TimeInterval = 1, - TimeStepStart = 2, - NUPCOL = 3 -}; + Balance(); + explicit Balance(const DeckKeyword& keyword); + explicit Balance(bool network_active); - Balance() = default; - Balance(bool network_active, const Tuning& tuning); - Balance(const Tuning& tuning, const DeckKeyword& keyword); + static Balance serializeObject(); CalcMode mode() const; double interval() const; @@ -51,11 +52,10 @@ enum class CalcMode { std::size_t pressure_max_iter() const; double thp_tolerance() const; std::size_t thp_max_iter() const; - std::optional target_balance_error() const; - std::optional max_balance_error() const; - double min_tstep() const; + const std::optional& target_balance_error() const; + const std::optional& max_balance_error() const; + const std::optional& min_tstep() const; - static Balance serializeObject(); bool operator==(const Balance& other) const; template @@ -72,7 +72,6 @@ enum class CalcMode { serializer(this->m_min_tstep); } - private: CalcMode calc_mode{CalcMode::None}; double calc_interval; @@ -82,11 +81,11 @@ private: double m_thp_tolerance; std::size_t m_thp_max_iter; - std::optional target_branch_balance_error; - std::optional max_branch_balance_error; - double m_min_tstep; + std::optional target_branch_balance_error{}; + std::optional max_branch_balance_error{}; + std::optional m_min_tstep{}; }; -} -} -#endif +}} // Opm::Network + +#endif // NETWORK_BALANCE_HPP diff --git a/opm/output/eclipse/DoubHEAD.hpp b/opm/output/eclipse/DoubHEAD.hpp index c4dd36ef3..b849c70bc 100644 --- a/opm/output/eclipse/DoubHEAD.hpp +++ b/opm/output/eclipse/DoubHEAD.hpp @@ -28,6 +28,7 @@ namespace Opm { struct Tuning; class Schedule; class UDQParams; + class UnitSystem; } namespace Opm { namespace RestartIO { @@ -57,7 +58,9 @@ namespace Opm { namespace RestartIO { double min_ec_grad; }; - struct NetBalanceParams { + struct NetBalanceParams { + explicit NetBalanceParams(const UnitSystem& usys); + double balancingInterval; double convTolNodPres; double convTolTHPCalc; diff --git a/opm/output/eclipse/VectorItems/doubhead.hpp b/opm/output/eclipse/VectorItems/doubhead.hpp index f7dea69c1..bed769a45 100644 --- a/opm/output/eclipse/VectorItems/doubhead.hpp +++ b/opm/output/eclipse/VectorItems/doubhead.hpp @@ -43,12 +43,21 @@ namespace Opm { namespace RestartIO { namespace Helpers { namespace VectorItems XxxMBE = 18, XxxLCV = 19, XxxWFL = 20, - Netbalint = 51, // balancingInterval - Netbalnpre = 53, // convTolNodPres - Netbalthpc = 50, // convTolTHPCalc - Netbaltarerr = 63, // targBranchBalError - Netbalmaxerr = 64, // maxBranchBalError - Netbalstepsz = 66, // minTimeStepSize + + Netbalthpc = 50, // Network balancing THP convergence limit (NETBALAN(4)) + Netbalint = 51, // Network balancing interval (NETBALAN(1)) + Netbalnpre = 53, // Network balancing nodal pressure + // convergence limit (NETBALAN(2)) + + Netbaltarerr = 63, // Target largest branch network balancing + // error at end of timestep (NETBALAN(6)) + + Netbalmaxerr = 64, // Maximum permitted network balancing error + // at end of timestep (NETBALAN(7)) + + Netbalstepsz = 66, // Minimum stepsize for steps limited by + // network balancing errors (NETBALAN(8)) + TrgDPR = 82, TfDiff = 83, DdpLim = 84, @@ -75,6 +84,14 @@ namespace Opm { namespace RestartIO { namespace Helpers { namespace VectorItems UdqPar_4 = 214, // UDQPARAM item number 4 (fractional equality tolerance used in ==, <= etc. functions) }; + namespace DoubHeadValue { + // Default if no active network (BRANPROP/NODEPROP) + constexpr auto NetBalNodPressDefault = 0.0; // Barsa + + // Default => Use TSMINZ from TUNING + constexpr auto NetBalMinTSDefault = 0.0; + } + }}}} // Opm::RestartIO::Helpers::VectorItems #endif // OPM_OUTPUT_ECLIPSE_VECTOR_DOUBHEAD_HPP diff --git a/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp b/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp index 64fda726e..60fa99aac 100644 --- a/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp +++ b/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp @@ -738,8 +738,8 @@ File {} line {}.)", wname, location.keyword, location.filename, location.lineno) } void Schedule::handleNETBALAN(HandlerContext& handlerContext) { - Network::Balance new_balance(this->snapshots.back().tuning(), handlerContext.keyword); - this->snapshots.back().network_balance.update( std::move(new_balance) ); + this->snapshots.back().network_balance + .update(Network::Balance{ handlerContext.keyword }); } void Schedule::handleNEXTSTEP(HandlerContext& handlerContext) { diff --git a/src/opm/input/eclipse/Schedule/Network/Balance.cpp b/src/opm/input/eclipse/Schedule/Network/Balance.cpp index 384ad2a0c..c97a9fca1 100644 --- a/src/opm/input/eclipse/Schedule/Network/Balance.cpp +++ b/src/opm/input/eclipse/Schedule/Network/Balance.cpp @@ -17,69 +17,103 @@ along with OPM. If not, see . */ - -#include #include -#include -#include + #include -namespace Opm { -namespace Network { +#include + +#include + +namespace { + Opm::Network::Balance::CalcMode + getNetworkBalancingMode(const double interval) + { + using CalcMode = Opm::Network::Balance::CalcMode; + + if (interval < 0.0) { + return CalcMode::NUPCOL; + } + else if (interval > 0.0) { + return CalcMode::TimeInterval; + } + else { + return CalcMode::TimeStepStart; + } + } + + double defaultCalcInterval() + { + static const auto interval = Opm::UnitSystem::newMETRIC() + .to_si(Opm::UnitSystem::measure::time, + Opm::ParserKeywords::NETBALAN::TIME_INTERVAL::defaultValue); + + return interval; + } + + double defaultNodePressureTolerance() + { + static const auto tol = Opm::UnitSystem::newMETRIC() + .to_si(Opm::UnitSystem::measure::pressure, + Opm::ParserKeywords::NETBALAN::PRESSURE_CONVERGENCE_LIMIT::defaultValue); + + return tol; + } +} + +namespace Opm { namespace Network { using NB = ParserKeywords::NETBALAN; -Balance::Balance(const Tuning& tuning, const DeckKeyword& keyword) { +Balance::Balance() + : calc_mode (CalcMode::None) + , calc_interval (defaultCalcInterval()) + , ptol (defaultNodePressureTolerance()) + , m_pressure_max_iter(NB::MAX_ITER::defaultValue) + , m_thp_tolerance (NB::THP_CONVERGENCE_LIMIT::defaultValue) + , m_thp_max_iter (NB::MAX_ITER_THP::defaultValue) +{} + +Balance::Balance(const DeckKeyword& keyword) +{ const auto& record = keyword[0]; - double interval = record.getItem().getSIDouble(0); - - if (interval < 0) { - this->calc_mode = CalcMode::NUPCOL; - this->calc_interval = 0.; - } - else if (interval == 0) { - this->calc_mode = CalcMode::TimeStepStart; - this->calc_interval = interval; - } - else { - this->calc_mode = CalcMode::TimeInterval; - this->calc_interval = interval; - } + this->calc_interval = record.getItem().getSIDouble(0); + this->calc_mode = getNetworkBalancingMode(this->calc_interval); this->ptol = record.getItem().getSIDouble(0); this->m_pressure_max_iter = record.getItem().get(0); this->m_thp_tolerance = record.getItem().get(0); this->m_thp_max_iter = record.getItem().get(0); - this->target_branch_balance_error = record.getItem().getSIDouble(0); - this->max_branch_balance_error = record.getItem().getSIDouble(0); - auto tstep_item = record.getItem(); - if (tstep_item.defaultApplied(0)) - this->m_min_tstep = tuning.TSMINZ; - else - this->m_min_tstep = record.getItem().getSIDouble(0); + if (const auto& targBE = record.getItem(); !targBE.defaultApplied(0)) { + this->target_branch_balance_error = targBE.getSIDouble(0); + } + + if (const auto& maxBE = record.getItem(); !maxBE.defaultApplied(0)) { + this->max_branch_balance_error = maxBE.getSIDouble(0); + } + + if (const auto& minTStep = record.getItem(); !minTStep.defaultApplied(0)) { + this->m_min_tstep = minTStep.getSIDouble(0); + } } - -Balance::Balance(bool network_active, const Tuning& tuning) - : calc_mode(CalcMode::TimeStepStart) - , calc_interval(0) +Balance::Balance(const bool network_active) + : calc_mode (CalcMode::TimeStepStart) + , calc_interval (defaultCalcInterval()) , m_thp_tolerance(NB::THP_CONVERGENCE_LIMIT::defaultValue) - , m_thp_max_iter( NB::MAX_ITER_THP::defaultValue ) + , m_thp_max_iter (NB::MAX_ITER_THP::defaultValue) { - NB parser_keyword{}; - UnitSystem default_units(UnitSystem::UnitType::UNIT_TYPE_METRIC); - if (network_active) { - const auto& item = parser_keyword.getRecord(0).get(NB::PRESSURE_CONVERGENCE_LIMIT::itemName); - this->ptol = default_units.to_si( item.dimensions()[0], item.getDefault()); + this->ptol = UnitSystem::newMETRIC().to_si(UnitSystem::measure::pressure, + NB::PRESSURE_CONVERGENCE_LIMIT::defaultValue); + this->m_pressure_max_iter = NB::MAX_ITER::defaultValue; - this->m_min_tstep = tuning.TSMINZ; - } else { - this->ptol = 0; + } + else { + this->ptol = 0.0; this->m_pressure_max_iter = 0; - this->m_min_tstep = 0; } } @@ -107,15 +141,17 @@ std::size_t Balance::pressure_max_iter() const { return this->m_pressure_max_iter; } -std::optional Balance::target_balance_error() const { +const std::optional& +Balance::target_balance_error() const { return this->target_branch_balance_error; } -std::optional Balance::max_balance_error() const { +const std::optional& +Balance::max_balance_error() const { return this->max_branch_balance_error; } -double Balance::min_tstep() const { +const std::optional& Balance::min_tstep() const { return this->m_min_tstep; } @@ -140,5 +176,4 @@ bool Balance::operator==(const Balance& other) const { this->m_min_tstep == other.m_min_tstep; } -} -} +}} // Opm::Network diff --git a/src/opm/input/eclipse/Schedule/Schedule.cpp b/src/opm/input/eclipse/Schedule/Schedule.cpp index e71a6134b..de0e6f7d5 100644 --- a/src/opm/input/eclipse/Schedule/Schedule.cpp +++ b/src/opm/input/eclipse/Schedule/Schedule.cpp @@ -2051,7 +2051,7 @@ void Schedule::create_first(const time_point& start_time, const std::optionalm_static.rst_config ) ); - sched_state.network_balance.update( Network::Balance(runspec.networkDimensions().active(), sched_state.tuning()) ); + sched_state.network_balance.update(Network::Balance{ runspec.networkDimensions().active() }); sched_state.update_sumthin(this->m_static.sumthin); sched_state.rptonly(this->m_static.rptonly); //sched_state.update_date( start_time ); diff --git a/src/opm/output/eclipse/CreateDoubHead.cpp b/src/opm/output/eclipse/CreateDoubHead.cpp index 09e6bf1da..9f8763d71 100644 --- a/src/opm/output/eclipse/CreateDoubHead.cpp +++ b/src/opm/output/eclipse/CreateDoubHead.cpp @@ -23,15 +23,19 @@ #include #include -#include + #include -#include +#include +#include + #include #include + #include #include #include +#include #include namespace { @@ -69,11 +73,10 @@ namespace { return static_cast(Opm::Metric::Time); } - - + Opm::RestartIO::DoubHEAD::guideRate computeGuideRate(const ::Opm::Schedule& sched, - const std::size_t lookup_step) + const std::size_t lookup_step) { double a = 0.; double b = 0.; @@ -112,7 +115,7 @@ namespace { Opm::RestartIO::DoubHEAD::liftOptPar computeLiftOptParam(const ::Opm::Schedule& sched, const Opm::UnitSystem& units, - const std::size_t lookup_step) + const std::size_t lookup_step) { using M = ::Opm::UnitSystem::measure; const auto& glo = sched.glo(lookup_step); @@ -122,40 +125,43 @@ namespace { }; } - Opm::RestartIO::DoubHEAD::NetBalanceParams - getNetworkBalanceParameters(const Opm::Schedule& sched, - const Opm::UnitSystem& units, - const std::size_t report_step) + void assignNetworkBalanceParameters(const Opm::Network::Balance& netbalan, + const Opm::UnitSystem& units, + Opm::RestartIO::DoubHEAD::NetBalanceParams& param) { using M = ::Opm::UnitSystem::measure; - double balancingInterval = 0.; - double convTolNodPres = 0.; - double convTolTHPCalc = 0.01; - double targBranchBalError = 1.E+20; - double maxBranchBalError = 1.E+20; - double minTimeStepSize = 0.; - if (report_step > 0) { - const auto& sched_state = sched[report_step]; - if (sched_state.network().active()) { - const auto lookup_step = report_step - 1; - balancingInterval = units.from_si(M::time, sched[lookup_step].network_balance().interval()); - convTolNodPres = units.from_si(M::pressure, sched[lookup_step].network_balance().pressure_tolerance()); - convTolTHPCalc = sched[lookup_step].network_balance().thp_tolerance(); - targBranchBalError = units.from_si(M::pressure, sched[lookup_step].network_balance().target_balance_error().value_or(Opm::ParserKeywords::NETBALAN::TARGET_BALANCE_ERROR::defaultValue)); - maxBranchBalError = units.from_si(M::pressure, sched[lookup_step].network_balance().max_balance_error().value_or(Opm::ParserKeywords::NETBALAN::MAX_BALANCE_ERROR::defaultValue)); - minTimeStepSize = units.from_si(M::time, sched[lookup_step].network_balance().min_tstep()); - } + param.balancingInterval = units.from_si(M::time, netbalan.interval()); + param.convTolNodPres = units.from_si(M::pressure, netbalan.pressure_tolerance()); + param.convTolTHPCalc = units.from_si(M::identity, netbalan.thp_tolerance()); + + if (const auto& trgBE = netbalan.target_balance_error(); trgBE.has_value()) { + param.targBranchBalError = units.from_si(M::pressure, trgBE.value()); } - return { - balancingInterval, - convTolNodPres, - convTolTHPCalc, - targBranchBalError, - maxBranchBalError, - minTimeStepSize - }; + if (const auto& maxBE = netbalan.max_balance_error(); maxBE.has_value()) { + param.maxBranchBalError = units.from_si(M::pressure, maxBE.value()); + } + + if (const auto& minTStep = netbalan.min_tstep(); minTStep.has_value()) { + param.minTimeStepSize = units.from_si(M::time, minTStep.value()); + } + } + + Opm::RestartIO::DoubHEAD::NetBalanceParams + getNetworkBalanceParameters(const Opm::Schedule& sched, + const Opm::UnitSystem& units, + const std::size_t report_step, + const std::size_t lookup_step) + { + auto param = Opm::RestartIO::DoubHEAD::NetBalanceParams{units}; + + if ((report_step > 0) && sched[lookup_step].network().active()) { + assignNetworkBalanceParameters(sched[lookup_step].network_balance(), + units, param); + } + + return param; } } // Anonymous @@ -173,17 +179,17 @@ createDoubHead(const EclipseState& es, const double nextTimeStep) { const auto& usys = es.getDeckUnitSystem(); - const auto& rspec = es.runspec(); + const auto& rspec = es.runspec(); const auto tconv = getTimeConv(usys); auto dh = DoubHEAD{} .tuningParameters(sched[lookup_step].tuning(), tconv) .timeStamp (computeTimeStamp(sched, simTime)) .drsdt (sched, lookup_step, tconv) - .udq_param(rspec.udqParams()) + .udq_param (rspec.udqParams()) .guide_rate_param(computeGuideRate(sched, lookup_step)) - .lift_opt_param(computeLiftOptParam(sched, usys, lookup_step)) - .netBalParams(getNetworkBalanceParameters(sched, usys, report_step)) + .lift_opt_param (computeLiftOptParam(sched, usys, lookup_step)) + .netBalParams (getNetworkBalanceParameters(sched, usys, report_step, lookup_step)) ; if (nextTimeStep > 0.0) { diff --git a/src/opm/output/eclipse/DoubHEAD.cpp b/src/opm/output/eclipse/DoubHEAD.cpp index 824699c1b..623dc0e67 100644 --- a/src/opm/output/eclipse/DoubHEAD.cpp +++ b/src/opm/output/eclipse/DoubHEAD.cpp @@ -18,20 +18,25 @@ */ #include -#include -#include -#include #include // Opm::RestartIO::makeUTCTime() - #include +#include +#include + +#include +#include + +#include + #include #include #include #include #include #include +#include #include #include #include @@ -368,12 +373,43 @@ namespace { // Failed to convert tm1 to time_t (unexpected). Use initial value. return toDateNum(tm0.tm_year, tm0.tm_yday); } + + double defaultNetBalanInterval() + { + static const auto interval = Opm::UnitSystem::newMETRIC() + .to_si(Opm::UnitSystem::measure::time, + Opm::ParserKeywords::NETBALAN::TIME_INTERVAL::defaultValue); + + return interval; + } + + double defaultNetBalanNodePressTol() + { + static const auto tol = Opm::UnitSystem::newMETRIC() + .to_si(Opm::UnitSystem::measure::pressure, + Opm::RestartIO::Helpers::VectorItems::DoubHeadValue::NetBalNodPressDefault); + + return tol; + } } // ===================================================================== // Public Interface (DoubHEAD member functions) Below Separator // --------------------------------------------------------------------- +Opm::RestartIO::DoubHEAD::NetBalanceParams::NetBalanceParams(const UnitSystem& usys) + : balancingInterval {usys.from_si(UnitSystem::measure::time, + defaultNetBalanInterval())} + , convTolNodPres {usys.from_si(UnitSystem::measure::pressure, + defaultNetBalanNodePressTol())} + , convTolTHPCalc {ParserKeywords::NETBALAN::THP_CONVERGENCE_LIMIT::defaultValue} // No usys.to_si() here! + , targBranchBalError{ParserKeywords::NETBALAN::TARGET_BALANCE_ERROR::defaultValue} // No usys.to_si() here! + , maxBranchBalError {ParserKeywords::NETBALAN::MAX_BALANCE_ERROR::defaultValue} // No usys.to_si() here! + , minTimeStepSize {Helpers::VectorItems::DoubHeadValue::NetBalMinTSDefault} +{} + +// --------------------------------------------------------------------- + Opm::RestartIO::DoubHEAD::DoubHEAD() : data_(Index::NUMBER_OF_ITEMS, 0.0) { @@ -650,13 +686,12 @@ Opm::RestartIO::DoubHEAD::lift_opt_param(const liftOptPar& lo_par) Opm::RestartIO::DoubHEAD& Opm::RestartIO::DoubHEAD::netBalParams(const NetBalanceParams& net_bal_par) { - this->data_[Netbalan_int] = net_bal_par.balancingInterval; + this->data_[Netbalan_int] = net_bal_par.balancingInterval; this->data_[Netbalan_npre] = net_bal_par.convTolNodPres; this->data_[Netbalan_thpc] = net_bal_par.convTolTHPCalc; - this->data_[Netbalan_tarerr] = net_bal_par.targBranchBalError; - this->data_[Netbalan_maxerr] = net_bal_par.maxBranchBalError; - this->data_[Netbalan_stepsz] = net_bal_par.minTimeStepSize; + this->data_[Netbalan_tarerr] = net_bal_par.targBranchBalError; + this->data_[Netbalan_maxerr] = net_bal_par.maxBranchBalError; + this->data_[Netbalan_stepsz] = net_bal_par.minTimeStepSize; return *this; } -