diff --git a/opm/simulators/flow/BlackoilModelEbosNldd.hpp b/opm/simulators/flow/BlackoilModelEbosNldd.hpp index 0503d8cb4..558ca23af 100644 --- a/opm/simulators/flow/BlackoilModelEbosNldd.hpp +++ b/opm/simulators/flow/BlackoilModelEbosNldd.hpp @@ -30,6 +30,7 @@ #include +#include #include #include #include @@ -331,7 +332,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(); @@ -695,26 +696,35 @@ private: const auto& solution = ebosSimulator.model().solution(0); std::vector domain_order(domains_.size()); - switch (model_.param().local_solve_approach_) { - case DomainSolveApproach::GaussSeidel: { + std::iota(domain_order.begin(), domain_order.end(), 0); + + 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()); 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] = avgpress; } - // 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. + 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] = maxpress; } break; } @@ -722,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) { @@ -730,28 +739,20 @@ 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] = maxres; } + break; } - break; - } - break; - } + } // end of switch (model_.param().local_domain_ordering_) - case DomainSolveApproach::Jacobi: - default: - std::iota(domain_order.begin(), domain_order.end(), 0); - break; + // Sort by largest measure, keeping index order if equal. + 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 @@ -764,7 +765,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 +789,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 @@ -883,8 +884,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); } diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index e6785ca9e..e40650376 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 { @@ -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}; @@ -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,15 +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 == "pressure") { - 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); } @@ -701,7 +693,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..27dc0dd35 100644 --- a/opm/simulators/flow/SubDomain.hpp +++ b/opm/simulators/flow/SubDomain.hpp @@ -22,6 +22,8 @@ #include +#include + #include #include @@ -36,9 +38,23 @@ namespace Opm //! \brief Measure to use for domain ordering. enum class DomainOrderingMeasure { AveragePressure, + MaxPressure, 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