OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #3032 from bska/netbalan-rst-fixes

Browse Source 
Increase Compatibility of NETBALAN Restart Output
This commit is contained in:
Markus Blatt 2022-05-27 21:46:21 +02:00 committed by GitHub
commit d278880091
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 228 additions and 133 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
opm/input/eclipse/Schedule/Network/Balance.hpp
View File

@ -17,33 +17,34 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef NETWORK_BALANCE_HPP
#define NETWORK_BALANCE_HPP
#include <optional>
#include <cstddef>
#include <cstddef>
#include <optional>
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<double> target_balance_error() const;
std::optional<double> max_balance_error() const;
double min_tstep() const;
const std::optional<double>& target_balance_error() const;
const std::optional<double>& max_balance_error() const;
const std::optional<double>& min_tstep() const;
static Balance serializeObject();
bool operator==(const Balance& other) const;
template<class Serializer>
@ -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<double> target_branch_balance_error;
std::optional<double> max_branch_balance_error;
double m_min_tstep;
std::optional<double> target_branch_balance_error{};
std::optional<double> max_branch_balance_error{};
std::optional<double> m_min_tstep{};
};
}
}
#endif
}} // Opm::Network
#endif // NETWORK_BALANCE_HPP

5
opm/output/eclipse/DoubHEAD.hpp
View File

@ -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;

29
opm/output/eclipse/VectorItems/doubhead.hpp
View File

@ -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

4
src/opm/input/eclipse/Schedule/KeywordHandlers.cpp
View File

@ -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) {

129
src/opm/input/eclipse/Schedule/Network/Balance.cpp
View File

@ -17,69 +17,103 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/input/eclipse/Parser/ParserKeywords/N.hpp>
#include <opm/input/eclipse/Schedule/Network/Balance.hpp>
#include <opm/input/eclipse/Schedule/Tuning.hpp>
#include <opm/input/eclipse/Deck/DeckKeyword.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
namespace Opm {
namespace Network {
#include <opm/input/eclipse/Parser/ParserKeywords/N.hpp>
#include <opm/input/eclipse/Deck/DeckKeyword.hpp>
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<NB::TIME_INTERVAL>().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<NB::TIME_INTERVAL>().getSIDouble(0);
this->calc_mode = getNetworkBalancingMode(this->calc_interval);
this->ptol = record.getItem<NB::PRESSURE_CONVERGENCE_LIMIT>().getSIDouble(0);
this->m_pressure_max_iter = record.getItem<NB::MAX_ITER>().get<int>(0);
this->m_thp_tolerance = record.getItem<NB::THP_CONVERGENCE_LIMIT>().get<double>(0);
this->m_thp_max_iter = record.getItem<NB::MAX_ITER_THP>().get<int>(0);
this->target_branch_balance_error = record.getItem<NB::TARGET_BALANCE_ERROR>().getSIDouble(0);
this->max_branch_balance_error = record.getItem<NB::MAX_BALANCE_ERROR>().getSIDouble(0);
auto tstep_item = record.getItem<NB::MIN_TIME_STEP>();
if (tstep_item.defaultApplied(0))
this->m_min_tstep = tuning.TSMINZ;
else
this->m_min_tstep = record.getItem<NB::MIN_TIME_STEP>().getSIDouble(0);
if (const auto& targBE = record.getItem<NB::TARGET_BALANCE_ERROR>(); !targBE.defaultApplied(0)) {
this->target_branch_balance_error = targBE.getSIDouble(0);
}
if (const auto& maxBE = record.getItem<NB::MAX_BALANCE_ERROR>(); !maxBE.defaultApplied(0)) {
this->max_branch_balance_error = maxBE.getSIDouble(0);
}
if (const auto& minTStep = record.getItem<NB::MIN_TIME_STEP>(); !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<double>());
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<double> Balance::target_balance_error() const {
const std::optional<double>&
Balance::target_balance_error() const {
return this->target_branch_balance_error;
}
std::optional<double> Balance::max_balance_error() const {
const std::optional<double>&
Balance::max_balance_error() const {
return this->max_branch_balance_error;
}
double Balance::min_tstep() const {
const std::optional<double>& 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

2
src/opm/input/eclipse/Schedule/Schedule.cpp
View File

@ -2051,7 +2051,7 @@ void Schedule::create_first(const time_point& start_time, const std::optional<ti
sched_state.guide_rate.update( GuideRateConfig() );
sched_state.rft_config.update( RFTConfig() );
sched_state.rst_config.update( RSTConfig::first( this->m_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 );

84
src/opm/output/eclipse/CreateDoubHead.cpp
View File

@ -23,15 +23,19 @@
#include <opm/output/eclipse/DoubHEAD.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Schedule/Group/GuideRateConfig.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/N.hpp>
#include <opm/input/eclipse/Schedule/Network/Balance.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/common/utility/TimeService.hpp>
#include <chrono>
#include <cstddef>
#include <optional>
#include <vector>
namespace {
@ -69,11 +73,10 @@ namespace {
return static_cast<double>(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) {

53
src/opm/output/eclipse/DoubHEAD.cpp
View File

@ -18,20 +18,25 @@
*/
#include <opm/output/eclipse/DoubHEAD.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Schedule/Tuning.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/output/eclipse/InteHEAD.hpp> // Opm::RestartIO::makeUTCTime()
#include <opm/output/eclipse/VectorItems/doubhead.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Schedule/Tuning.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/N.hpp>
#include <chrono>
#include <cmath>
#include <ctime>
#include <iterator>
#include <map>
#include <numeric>
#include <optional>
#include <ratio>
#include <utility>
#include <vector>
@ -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;
}