Parameterize network sub-iterations and pressure update dampening

This commit is contained in:
Vegard Kippe 2025-01-27 11:44:10 +01:00
parent 6bfb60ded0
commit 1d2e413d63
5 changed files with 33 additions and 10 deletions

View File

@ -97,6 +97,9 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
deck_file_name_ = Parameters::Get<Parameters::EclDeckFileName>();
network_max_strict_iterations_ = Parameters::Get<Parameters::NetworkMaxStrictIterations>();
network_max_iterations_ = Parameters::Get<Parameters::NetworkMaxIterations>();
network_max_sub_iterations_ = Parameters::Get<Parameters::NetworkMaxSubIterations>();
network_pressure_update_damping_factor_ = Parameters::Get<Parameters::NetworkPressureUpdateDampingFactor<Scalar>>();
network_max_pressure_update_in_bars_ = Parameters::Get<Parameters::NetworkMaxPressureUpdateInBars<Scalar>>();
local_domain_ordering_ = domainOrderingMeasureFromString(Parameters::Get<Parameters::LocalDomainsOrderingMeasure>());
write_partitions_ = Parameters::Get<Parameters::DebugEmitCellPartition>();
@ -220,6 +223,12 @@ void BlackoilModelParameters<Scalar>::registerParameters()
("Maximum iterations in network solver before relaxing tolerance");
Parameters::Register<Parameters::NetworkMaxIterations>
("Maximum number of iterations in the network solver before giving up");
Parameters::Register<Parameters::NetworkMaxSubIterations>
("Maximum number of sub-iterations to update network pressures (within a single well/group control update)");
Parameters::Register<Parameters::NetworkPressureUpdateDampingFactor<Scalar>>
("Damping factor in the inner network pressure update iterations");
Parameters::Register<Parameters::NetworkMaxPressureUpdateInBars<Scalar>>
("Maximum pressure update in the inner network pressure update iterations");
Parameters::Register<Parameters::NonlinearSolver>
("Choose nonlinear solver. Valid choices are newton or nldd.");
Parameters::Register<Parameters::LocalSolveApproach>

View File

@ -130,6 +130,11 @@ struct CheckGroupConstraintsInnerWellIterations { static constexpr bool value =
// Network solver parameters
struct NetworkMaxStrictIterations { static constexpr int value = 10; };
struct NetworkMaxIterations { static constexpr int value = 20; };
struct NetworkMaxSubIterations { static constexpr int value = 20; };
template<class Scalar>
struct NetworkPressureUpdateDampingFactor { static constexpr Scalar value = 0.1; };
template<class Scalar>
struct NetworkMaxPressureUpdateInBars { static constexpr Scalar value = 5.0; };
struct NonlinearSolver { static constexpr auto value = "newton"; };
struct LocalSolveApproach { static constexpr auto value = "gauss-seidel"; };
struct MaxLocalSolveIterations { static constexpr int value = 20; };
@ -297,6 +302,15 @@ public:
/// Maximum number of iterations in the network solver before giving up
int network_max_iterations_;
/// Maximum number of sub-iterations to update network pressures (within a single well/group control update)
int network_max_sub_iterations_;
/// Damping factor in the inner network pressure update iterations
Scalar network_pressure_update_damping_factor_;
/// Maximum pressure update in the inner network pressure update iterations
Scalar network_max_pressure_update_in_bars_;
/// Nonlinear solver type: newton or nldd.
std::string nonlinear_solver_;
/// 'jacobi' and 'gauss-seidel' supported.

View File

@ -1516,7 +1516,7 @@ inferLocalShutWells()
template<class Scalar>
Scalar BlackoilWellModelGeneric<Scalar>::
updateNetworkPressures(const int reportStepIdx, const Scalar damping_factor)
updateNetworkPressures(const int reportStepIdx, const Scalar damping_factor, const Scalar upper_update_bound)
{
OPM_TIMEFUNCTION();
// Get the network and return if inactive (no wells in network at this time)
@ -1552,11 +1552,9 @@ updateNetworkPressures(const int reportStepIdx, const Scalar damping_factor)
if (std::abs(change) > network_imbalance) {
network_imbalance = std::abs(change);
}
// we dampen the amount of the nodal pressure can change during one iteration
// due to the fact our nodal pressure calculation is somewhat explicit
// TODO: the following parameters are subject to adjustment for optimization purpose
constexpr Scalar upper_update_bound = 5.0 * unit::barsa;
// relative dampening factor based on update value
// 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;

View File

@ -351,7 +351,8 @@ protected:
bool wasDynamicallyShutThisTimeStep(const int well_index) const;
Scalar updateNetworkPressures(const int reportStepIdx,
const Scalar damping_factor);
const Scalar damping_factor,
const Scalar update_upper_bound);
void updateWsolvent(const Group& group,
const int reportStepIdx,

View File

@ -1919,11 +1919,12 @@ namespace Opm {
const double dt = this->simulator_.timeStepSize();
// Calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES)
bool well_group_thp_updated = computeWellGroupThp(dt, deferred_logger);
constexpr int max_number_of_sub_iterations = 20;
constexpr Scalar damping_factor = 0.1;
const int max_number_of_sub_iterations = param_.network_max_sub_iterations_;
const Scalar network_pressure_update_damping_factor = param_.network_pressure_update_damping_factor_;
const Scalar network_max_pressure_update = param_.network_max_pressure_update_in_bars_ * unit::barsa;
bool more_network_sub_update = false;
for (int i = 0; i < max_number_of_sub_iterations; i++) {
const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx, damping_factor);
const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx, network_pressure_update_damping_factor, network_max_pressure_update);
const Scalar network_imbalance = comm.max(local_network_imbalance);
const auto& balance = this->schedule()[episodeIdx].network_balance();
constexpr Scalar relaxation_factor = 10.0;