Merge pull request #5021 from atgeirr/nldd-refinements

Nldd refinements

opm/simulators/flow/BlackoilModelEbosNldd.hpp

@@ -30,6 +30,7 @@
 
 #include <opm/simulators/aquifers/AquiferGridUtils.hpp>
 
+#include <opm/simulators/flow/countGlobalCells.hpp>
 #include <opm/simulators/flow/partitionCells.hpp>
 #include <opm/simulators/flow/priVarsPacking.hpp>
 #include <opm/simulators/flow/SubDomain.hpp>
@@ -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<int> domain_order(domains_.size());
-        switch (model_.param().local_solve_approach_) {
+        std::iota(domain_order.begin(), domain_order.end(), 0);
-        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<double> measure_per_domain(domains_.size());
             switch (model_.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);
+                    measure_per_domain[domain.index] = avgpress;
                 }
-                // Lexicographical sort by pressure, then index.
+                break;
-                std::sort(avgpress_per_domain.begin(), avgpress_per_domain.end());
+            }
-                // Reverse since we want high-pressure regions solved first.
+            case DomainOrderingMeasure::MaxPressure: {
-                std::reverse(avgpress_per_domain.begin(), avgpress_per_domain.end());
+                // Use max pressures to order domains.
-                for (std::size_t ii = 0; ii < domains_.size(); ++ii) {
+                for (const auto& domain : domains_) {
-                    domain_order[ii] = avgpress_per_domain[ii].second;
+                    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<std::pair<double, int>> 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);
+                    measure_per_domain[domain.index] = maxres;
-                }
-                // 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;
                 }
+                break;
             }
-            break;
+            } // end of switch (model_.param().local_domain_ordering_)
-            }
-            break;
-        }
 
-        case DomainSolveApproach::Jacobi:
+            // Sort by largest measure, keeping index order if equal.
-        default:
+            const auto& m = measure_per_domain;
-            std::iota(domain_order.begin(), domain_order.end(), 0);
+            std::stable_sort(domain_order.begin(), domain_order.end(),
-            break;
+                             [&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<class GlobalEqVector>
@@ -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<Well>{};
 
+        // 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);
     }
 

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -402,7 +402,7 @@ struct NonlinearSolver<TypeTag, TTag::FlowModelParameters> {
 };
 template<class TypeTag>
 struct LocalSolveApproach<TypeTag, TTag::FlowModelParameters> {
-    static constexpr auto value = "jacobi";
+    static constexpr auto value = "gauss-seidel";
 };
 template<class TypeTag>
 struct MaxLocalSolveIterations<TypeTag, TTag::FlowModelParameters> {
@@ -416,7 +416,7 @@ struct LocalToleranceScalingMb<TypeTag, TTag::FlowModelParameters> {
 template<class TypeTag>
 struct LocalToleranceScalingCnv<TypeTag, TTag::FlowModelParameters> {
     using type = GetPropType<TypeTag, Scalar>;
-    static constexpr type value = 0.01;
+    static constexpr type value = 0.1;
 };
 template<class TypeTag>
 struct NlddNumInitialNewtonIter<TypeTag, TTag::FlowModelParameters> {
@@ -439,7 +439,7 @@ struct LocalDomainsPartitioningMethod<TypeTag, TTag::FlowModelParameters> {
 };
 template<class TypeTag>
 struct LocalDomainsOrderingMeasure<TypeTag, TTag::FlowModelParameters> {
-    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);
+            local_domain_ordering_ = domainOrderingMeasureFromString(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.");
-            }
-
             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.");
 

opm/simulators/flow/SubDomain.hpp

@@ -22,6 +22,8 @@
 
 #include <opm/grid/common/SubGridPart.hpp>
 
+#include <fmt/format.h>
+
 #include <utility>
 #include <vector>
 
@@ -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 <class Grid>