Merge pull request #4738 from akva2/blackoilmodelebos_domain_ordering

BlackoilModelEbos: improve domain ordering

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<int> domain_order(domains_.size());
+            const auto domain_order = this->getSubdomainOrder();
-            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<std::pair<double, int>> 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<std::pair<double, int>> 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);
-            }
 
             // -----------   Solve each domain separately   -----------
             std::vector<SimulatorReportSingle> 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<int> getSubdomainOrder()
+        {
+            const auto& solution = ebosSimulator().model().solution(0);
+
+            std::vector<int> 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<std::pair<double, int>> 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<std::pair<double, int>> 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<bool> wasSwitched_;
     };

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -20,9 +20,15 @@
 #ifndef OPM_BLACKOILMODELPARAMETERS_EBOS_HEADER_INCLUDED
 #define OPM_BLACKOILMODELPARAMETERS_EBOS_HEADER_INCLUDED
 
-#include <opm/models/utils/propertysystem.hh>
+#include <opm/models/discretization/common/fvbaseproperties.hh>
-#include <opm/models/utils/parametersystem.hh>
 
+#include <opm/models/utils/basicproperties.hh>
+#include <opm/models/utils/parametersystem.hh>
+#include <opm/models/utils/propertysystem.hh>
+
+#include <opm/simulators/flow/SubDomain.hpp>
+
+#include <stdexcept>
 #include <string>
 
 namespace Opm::Properties {
@@ -209,6 +215,10 @@ template<class TypeTag, class MyTypeTag>
 struct LocalDomainsPartitioningMethod {
     using type = UndefinedProperty;
 };
+template<class TypeTag, class MyTypeTag>
+struct LocalDomainsOrderingMeasure {
+    using type = UndefinedProperty;
+};
 template<class TypeTag>
 struct DbhpMaxRel<TypeTag, TTag::FlowModelParameters> {
     using type = GetPropType<TypeTag, Scalar>;
@@ -401,6 +411,10 @@ template<class TypeTag>
 struct LocalDomainsPartitioningMethod<TypeTag, TTag::FlowModelParameters> {
     static constexpr auto value = "zoltan";
 };
+template<class TypeTag>
+struct LocalDomainsOrderingMeasure<TypeTag, TTag::FlowModelParameters> {
+    static constexpr auto value = "pressure";
+};
 // if openMP is available, determine the number threads per process automatically.
 #if _OPENMP
 template<class TypeTag>
@@ -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

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 <class Grid>