From 9c9fae26a82d9a8bae0b50235d8c3b48fb82029d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 15 Nov 2023 09:35:08 +0100 Subject: [PATCH 1/6] Require initial assembly for subdomains with Gauss-Seidel. --- opm/simulators/flow/BlackoilModelEbosNldd.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbosNldd.hpp b/opm/simulators/flow/BlackoilModelEbosNldd.hpp index 0503d8cb4..695b13935 100644 --- a/opm/simulators/flow/BlackoilModelEbosNldd.hpp +++ b/opm/simulators/flow/BlackoilModelEbosNldd.hpp @@ -331,7 +331,7 @@ private: solveDomain(const Domain& domain, const SimulatorTimerInterface& timer, [[maybe_unused]] const int global_iteration, - const bool initial_assembly_required = false) + const bool initial_assembly_required) { auto& ebosSimulator = model_.ebosSimulator(); @@ -764,7 +764,7 @@ private: { auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain); auto initial_local_solution = Details::extractVector(solution, domain.cells); - auto res = solveDomain(domain, timer, iteration); + auto res = solveDomain(domain, timer, iteration, false); local_report = res.first; if (local_report.converged) { auto local_solution = Details::extractVector(solution, domain.cells); @@ -788,7 +788,7 @@ private: { auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain); auto initial_local_solution = Details::extractVector(solution, domain.cells); - auto res = solveDomain(domain, timer, iteration); + auto res = solveDomain(domain, timer, iteration, true); local_report = res.first; if (!local_report.converged) { // We look at the detailed convergence report to evaluate From 8acc8af2a41ced9d3391255c52a2821bba80bc8c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 15 Nov 2023 09:36:01 +0100 Subject: [PATCH 2/6] Add MaxPressure NLDD domain ordering option, make it default. --- opm/simulators/flow/BlackoilModelEbosNldd.hpp | 47 ++++++++++++------- .../flow/BlackoilModelParametersEbos.hpp | 10 ++-- opm/simulators/flow/SubDomain.hpp | 1 + 3 files changed, 36 insertions(+), 22 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbosNldd.hpp b/opm/simulators/flow/BlackoilModelEbosNldd.hpp index 695b13935..607a89297 100644 --- a/opm/simulators/flow/BlackoilModelEbosNldd.hpp +++ b/opm/simulators/flow/BlackoilModelEbosNldd.hpp @@ -695,26 +695,33 @@ private: const auto& solution = ebosSimulator.model().solution(0); std::vector domain_order(domains_.size()); + switch (model_.param().local_solve_approach_) { case DomainSolveApproach::GaussSeidel: { + // Calculate the measure used to order the domains. + std::vector> measure_per_domain(domains_.size()); switch (model_.param().local_domain_ordering_) { case DomainOrderingMeasure::AveragePressure: { // Use average pressures to order domains. - std::vector> avgpress_per_domain(domains_.size()); for (const auto& domain : domains_) { double press_sum = 0.0; for (const int c : domain.cells) { press_sum += solution[c][Indices::pressureSwitchIdx]; } const double avgpress = press_sum / domain.cells.size(); - avgpress_per_domain[domain.index] = std::make_pair(avgpress, domain.index); + measure_per_domain[domain.index] = std::make_pair(avgpress, domain.index); } - // Lexicographical sort by pressure, then index. - std::sort(avgpress_per_domain.begin(), avgpress_per_domain.end()); - // Reverse since we want high-pressure regions solved first. - std::reverse(avgpress_per_domain.begin(), avgpress_per_domain.end()); - for (std::size_t ii = 0; ii < domains_.size(); ++ii) { - domain_order[ii] = avgpress_per_domain[ii].second; + break; + } + case DomainOrderingMeasure::MaxPressure: { + // Use max pressures to order domains. + std::vector> maxpress_per_domain(domains_.size()); + for (const auto& domain : domains_) { + double maxpress = 0.0; + for (const int c : domain.cells) { + maxpress = std::max(maxpress, solution[c][Indices::pressureSwitchIdx]); + } + measure_per_domain[domain.index] = std::make_pair(maxpress, domain.index); } break; } @@ -730,20 +737,24 @@ private: maxres = std::max(maxres, std::fabs(residual[c][ii])); } } - maxres_per_domain[domain.index] = std::make_pair(maxres, domain.index); - } - // Lexicographical sort by pressure, then index. - std::sort(maxres_per_domain.begin(), maxres_per_domain.end()); - // Reverse since we want high-pressure regions solved first. - std::reverse(maxres_per_domain.begin(), maxres_per_domain.end()); - for (std::size_t ii = 0; ii < domains_.size(); ++ii) { - domain_order[ii] = maxres_per_domain[ii].second; + measure_per_domain[domain.index] = std::make_pair(maxres, domain.index); } + break; } - break; + } // end of switch (model_.param().local_domain_ordering_) + + + // Sort by largest measure, keeping index order if equal. + std::stable_sort(measure_per_domain.begin(), measure_per_domain.end(), + [](const auto& m1, const auto& m2){ return m1.first > m2.first; }); + + // Assign domain order. + for (std::size_t ii = 0; ii < domains_.size(); ++ii) { + domain_order[ii] = measure_per_domain[ii].second; } + break; - } + } // end of case DomainSolveApproach::GaussSeidel case DomainSolveApproach::Jacobi: default: diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index e6785ca9e..4a6d9f136 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -439,7 +439,7 @@ struct LocalDomainsPartitioningMethod { }; template struct LocalDomainsOrderingMeasure { - static constexpr auto value = "pressure"; + static constexpr auto value = "maxpressure"; }; // if openMP is available, determine the number threads per process automatically. #if _OPENMP @@ -579,7 +579,7 @@ namespace Opm int num_local_domains_{0}; double local_domain_partition_imbalance_{1.03}; std::string local_domain_partition_method_; - DomainOrderingMeasure local_domain_ordering_{DomainOrderingMeasure::AveragePressure}; + DomainOrderingMeasure local_domain_ordering_{DomainOrderingMeasure::MaxPressure}; bool write_partitions_{false}; @@ -642,7 +642,9 @@ namespace Opm std::string measure = EWOMS_GET_PARAM(TypeTag, std::string, LocalDomainsOrderingMeasure); if (measure == "residual") { local_domain_ordering_ = DomainOrderingMeasure::Residual; - } else if (measure == "pressure") { + } else if (measure == "maxpressure") { + local_domain_ordering_ = DomainOrderingMeasure::MaxPressure; + } else if (measure == "averagepressure") { local_domain_ordering_ = DomainOrderingMeasure::AveragePressure; } else { throw std::runtime_error("Invalid domain ordering '" + measure + "' specified."); @@ -701,7 +703,7 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, std::string, LocalDomainsPartitioningMethod, "Subdomain partitioning method. " "Allowed values are 'zoltan', 'simple', and the name of a partition file ending with '.partition'."); EWOMS_REGISTER_PARAM(TypeTag, std::string, LocalDomainsOrderingMeasure, "Subdomain ordering measure. " - "Allowed values are 'pressure' and 'residual'."); + "Allowed values are 'maxpressure', 'averagepressure' and 'residual'."); EWOMS_REGISTER_PARAM(TypeTag, bool, DebugEmitCellPartition, "Whether or not to emit cell partitions as a debugging aid."); diff --git a/opm/simulators/flow/SubDomain.hpp b/opm/simulators/flow/SubDomain.hpp index 033d4f772..ba92c9ff8 100644 --- a/opm/simulators/flow/SubDomain.hpp +++ b/opm/simulators/flow/SubDomain.hpp @@ -36,6 +36,7 @@ namespace Opm //! \brief Measure to use for domain ordering. enum class DomainOrderingMeasure { AveragePressure, + MaxPressure, Residual }; From 94554ea4f243633b8ec944ad255523c4c9642e19 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Nov 2023 10:37:43 +0100 Subject: [PATCH 3/6] Change default nldd number of domains to one per 1000 cells. --- opm/simulators/flow/BlackoilModelEbosNldd.hpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/opm/simulators/flow/BlackoilModelEbosNldd.hpp b/opm/simulators/flow/BlackoilModelEbosNldd.hpp index 607a89297..046c322d1 100644 --- a/opm/simulators/flow/BlackoilModelEbosNldd.hpp +++ b/opm/simulators/flow/BlackoilModelEbosNldd.hpp @@ -30,6 +30,7 @@ #include +#include #include #include #include @@ -894,8 +895,15 @@ private: ? this->model_.ebosSimulator().vanguard().schedule().getWellsatEnd() : std::vector{}; + // If defaulted parameter for number of domains, choose a reasonable default. + constexpr int default_cells_per_domain = 1000; + const int num_cells = Opm::detail::countGlobalCells(grid); + const int num_domains = param.num_local_domains_ > 0 + ? param.num_local_domains_ + : num_cells / default_cells_per_domain; + return ::Opm::partitionCells(param.local_domain_partition_method_, - param.num_local_domains_, + num_domains, grid.leafGridView(), wells, zoltan_ctrl); } From 7578fbf1440e47866f3877a9a0fe1c91d6de26ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Nov 2023 10:38:11 +0100 Subject: [PATCH 4/6] Change nldd default approach to gauss-seidel. Also adjust local CNV tolerance scaling to 0.1 instead of 0.01. --- opm/simulators/flow/BlackoilModelParametersEbos.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index 4a6d9f136..25131f5fe 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -402,7 +402,7 @@ struct NonlinearSolver { }; template struct LocalSolveApproach { - static constexpr auto value = "jacobi"; + static constexpr auto value = "gauss-seidel"; }; template struct MaxLocalSolveIterations { @@ -416,7 +416,7 @@ struct LocalToleranceScalingMb { template struct LocalToleranceScalingCnv { using type = GetPropType; - static constexpr type value = 0.01; + static constexpr type value = 0.1; }; template struct NlddNumInitialNewtonIter { From 6f04c31c7c61f161c42839e407a6d8f32b9d67ea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Nov 2023 13:37:22 +0100 Subject: [PATCH 5/6] Use string to enum helper. --- .../flow/BlackoilModelParametersEbos.hpp | 14 ++------------ opm/simulators/flow/SubDomain.hpp | 15 +++++++++++++++ 2 files changed, 17 insertions(+), 12 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index 25131f5fe..e40650376 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -620,7 +620,7 @@ namespace Opm use_average_density_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseAverageDensityMsWells); local_well_solver_control_switching_ = EWOMS_GET_PARAM(TypeTag, bool, LocalWellSolveControlSwitching); nonlinear_solver_ = EWOMS_GET_PARAM(TypeTag, std::string, NonlinearSolver); - std::string approach = EWOMS_GET_PARAM(TypeTag, std::string, LocalSolveApproach); + const auto approach = EWOMS_GET_PARAM(TypeTag, std::string, LocalSolveApproach); if (approach == "jacobi") { local_solve_approach_ = DomainSolveApproach::Jacobi; } else if (approach == "gauss-seidel") { @@ -639,17 +639,7 @@ namespace Opm deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName); network_max_strict_iterations_ = EWOMS_GET_PARAM(TypeTag, int, NetworkMaxStrictIterations); network_max_iterations_ = EWOMS_GET_PARAM(TypeTag, int, NetworkMaxIterations); - std::string measure = EWOMS_GET_PARAM(TypeTag, std::string, LocalDomainsOrderingMeasure); - if (measure == "residual") { - local_domain_ordering_ = DomainOrderingMeasure::Residual; - } else if (measure == "maxpressure") { - local_domain_ordering_ = DomainOrderingMeasure::MaxPressure; - } else if (measure == "averagepressure") { - local_domain_ordering_ = DomainOrderingMeasure::AveragePressure; - } else { - throw std::runtime_error("Invalid domain ordering '" + measure + "' specified."); - } - + local_domain_ordering_ = domainOrderingMeasureFromString(EWOMS_GET_PARAM(TypeTag, std::string, LocalDomainsOrderingMeasure)); write_partitions_ = EWOMS_GET_PARAM(TypeTag, bool, DebugEmitCellPartition); } diff --git a/opm/simulators/flow/SubDomain.hpp b/opm/simulators/flow/SubDomain.hpp index ba92c9ff8..27dc0dd35 100644 --- a/opm/simulators/flow/SubDomain.hpp +++ b/opm/simulators/flow/SubDomain.hpp @@ -22,6 +22,8 @@ #include +#include + #include #include @@ -40,6 +42,19 @@ namespace Opm Residual }; + inline DomainOrderingMeasure domainOrderingMeasureFromString(const std::string_view measure) + { + if (measure == "residual") { + return DomainOrderingMeasure::Residual; + } else if (measure == "maxpressure") { + return DomainOrderingMeasure::MaxPressure; + } else if (measure == "averagepressure") { + return DomainOrderingMeasure::AveragePressure; + } else { + throw std::runtime_error(fmt::format("Invalid domain ordering '{}' specified", measure)); + } + } + /// Representing a part of a grid, in a way suitable for performing /// local solves. template From 6945b927ec3fd0adb1ad31408854b489ee9c1545 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 Nov 2023 14:57:45 +0100 Subject: [PATCH 6/6] Minor refactoring of domain ordering for NLDD. --- opm/simulators/flow/BlackoilModelEbosNldd.hpp | 41 +++++++------------ 1 file changed, 15 insertions(+), 26 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbosNldd.hpp b/opm/simulators/flow/BlackoilModelEbosNldd.hpp index 046c322d1..558ca23af 100644 --- a/opm/simulators/flow/BlackoilModelEbosNldd.hpp +++ b/opm/simulators/flow/BlackoilModelEbosNldd.hpp @@ -696,11 +696,14 @@ private: const auto& solution = ebosSimulator.model().solution(0); std::vector domain_order(domains_.size()); + std::iota(domain_order.begin(), domain_order.end(), 0); - switch (model_.param().local_solve_approach_) { - case DomainSolveApproach::GaussSeidel: { + if (model_.param().local_solve_approach_ == DomainSolveApproach::Jacobi) { + // Do nothing, 0..n-1 order is fine. + return domain_order; + } else if (model_.param().local_solve_approach_ == DomainSolveApproach::GaussSeidel) { // Calculate the measure used to order the domains. - std::vector> measure_per_domain(domains_.size()); + std::vector measure_per_domain(domains_.size()); switch (model_.param().local_domain_ordering_) { case DomainOrderingMeasure::AveragePressure: { // Use average pressures to order domains. @@ -710,19 +713,18 @@ private: press_sum += solution[c][Indices::pressureSwitchIdx]; } const double avgpress = press_sum / domain.cells.size(); - measure_per_domain[domain.index] = std::make_pair(avgpress, domain.index); + measure_per_domain[domain.index] = avgpress; } break; } case DomainOrderingMeasure::MaxPressure: { // Use max pressures to order domains. - std::vector> maxpress_per_domain(domains_.size()); for (const auto& domain : domains_) { double maxpress = 0.0; for (const int c : domain.cells) { maxpress = std::max(maxpress, solution[c][Indices::pressureSwitchIdx]); } - measure_per_domain[domain.index] = std::make_pair(maxpress, domain.index); + measure_per_domain[domain.index] = maxpress; } break; } @@ -730,7 +732,6 @@ private: // Use maximum residual to order domains. const auto& residual = ebosSimulator.model().linearizer().residual(); const int num_vars = residual[0].size(); - std::vector> maxres_per_domain(domains_.size()); for (const auto& domain : domains_) { double maxres = 0.0; for (const int c : domain.cells) { @@ -738,32 +739,20 @@ private: maxres = std::max(maxres, std::fabs(residual[c][ii])); } } - measure_per_domain[domain.index] = std::make_pair(maxres, domain.index); + measure_per_domain[domain.index] = maxres; } break; } } // end of switch (model_.param().local_domain_ordering_) - // Sort by largest measure, keeping index order if equal. - std::stable_sort(measure_per_domain.begin(), measure_per_domain.end(), - [](const auto& m1, const auto& m2){ return m1.first > m2.first; }); - - // Assign domain order. - for (std::size_t ii = 0; ii < domains_.size(); ++ii) { - domain_order[ii] = measure_per_domain[ii].second; - } - - break; - } // end of case DomainSolveApproach::GaussSeidel - - case DomainSolveApproach::Jacobi: - default: - std::iota(domain_order.begin(), domain_order.end(), 0); - break; + const auto& m = measure_per_domain; + std::stable_sort(domain_order.begin(), domain_order.end(), + [&m](const int i1, const int i2){ return m[i1] > m[i2]; }); + return domain_order; + } else { + throw std::logic_error("Domain solve approach must be Jacobi or Gauss-Seidel"); } - - return domain_order; } template