diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index b978f0a1e..e26183ab5 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -534,52 +534,7 @@ namespace Opm { auto locally_solved = initial_solution; // ----------- Decide on an ordering for the domains ----------- - std::vector domain_order(domains_.size()); - if (param_.local_solve_approach_ == "gauss-seidel") { - // TODO: enable flexibility and choice in choosing domain ordering approach. - if (true) { - // 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); - } - // 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 (size_t ii = 0; ii < domains_.size(); ++ii) { - domain_order[ii] = avgpress_per_domain[ii].second; - } - } else { - // 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) { - for (int ii = 0; ii < num_vars; ++ii) { - 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 (size_t ii = 0; ii < domains_.size(); ++ii) { - domain_order[ii] = maxres_per_domain[ii].second; - } - } - } else { - std::iota(domain_order.begin(), domain_order.end(), 0); - } + const auto domain_order = this->getSubdomainOrder(); // ----------- Solve each domain separately ----------- std::vector domain_reports(domains_.size()); @@ -1833,6 +1788,64 @@ namespace Opm { double maxResidualAllowed() const { return param_.max_residual_allowed_; } double linear_solve_setup_time_; + //! \brief Returns subdomain ordered according to method and ordering measure. + std::vector getSubdomainOrder() + { + const auto& solution = ebosSimulator().model().solution(0); + + std::vector domain_order(domains_.size()); + if (param_.local_solve_approach_ == "gauss-seidel") { + switch (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); + } + // 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 (size_t ii = 0; ii < domains_.size(); ++ii) { + domain_order[ii] = avgpress_per_domain[ii].second; + } + break; + } + case DomainOrderingMeasure::Residual: { + // 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) { + for (int ii = 0; ii < num_vars; ++ii) { + 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 (size_t ii = 0; ii < domains_.size(); ++ii) { + domain_order[ii] = maxres_per_domain[ii].second; + } + } + } + } else { + std::iota(domain_order.begin(), domain_order.end(), 0); + } + + return domain_order; + } + public: std::vector wasSwitched_; }; diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index f52b62e71..ba4e3f94d 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -20,9 +20,15 @@ #ifndef OPM_BLACKOILMODELPARAMETERS_EBOS_HEADER_INCLUDED #define OPM_BLACKOILMODELPARAMETERS_EBOS_HEADER_INCLUDED -#include -#include +#include +#include +#include +#include + +#include + +#include #include namespace Opm::Properties { @@ -209,6 +215,10 @@ template struct LocalDomainsPartitioningMethod { using type = UndefinedProperty; }; +template +struct LocalDomainsOrderingMeasure { + using type = UndefinedProperty; +}; template struct DbhpMaxRel { using type = GetPropType; @@ -401,6 +411,10 @@ template struct LocalDomainsPartitioningMethod { static constexpr auto value = "zoltan"; }; +template +struct LocalDomainsOrderingMeasure { + static constexpr auto value = "pressure"; +}; // if openMP is available, determine the number threads per process automatically. #if _OPENMP template @@ -535,6 +549,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}; /// Construct from user parameters or defaults. BlackoilModelParametersEbos() @@ -581,7 +596,15 @@ namespace Opm local_domain_partition_method_ = EWOMS_GET_PARAM(TypeTag, std::string, LocalDomainsPartitioningMethod); 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); + 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."); + } } static void registerParameters() @@ -631,6 +654,8 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalDomainsPartitioningImbalance, "Subdomain partitioning imbalance tolerance. 1.03 is 3 percent imbalance."); 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'."); } }; } // namespace Opm diff --git a/opm/simulators/flow/SubDomain.hpp b/opm/simulators/flow/SubDomain.hpp index 7796ff07e..e79297788 100644 --- a/opm/simulators/flow/SubDomain.hpp +++ b/opm/simulators/flow/SubDomain.hpp @@ -28,6 +28,12 @@ namespace Opm { + //! \brief Measure to use for domain ordering. + enum class DomainOrderingMeasure { + AveragePressure, + Residual + }; + /// Representing a part of a grid, in a way suitable for performing /// local solves. template