Merge pull request #3764 from atgeirr/aspin

Add nonlinear solver based on domain decomposition method

opm/simulators/flow/Banners.cpp

@@ -112,13 +112,20 @@ void printFlowBanner(int nprocs, int nthreads, std::string_view moduleVersionNam
 }
 
 void printFlowTrailer(int nprocs, int nthreads,
-                      const SimulatorReport& report)
+                      const SimulatorReport& report,
+                      const SimulatorReportSingle& localsolves_report)
 {
     std::ostringstream ss;
     ss << "\n\n================    End of simulation     ===============\n\n";
     ss << fmt::format("Number of MPI processes: {:9}\n", nprocs);
     ss << fmt::format("Threads per MPI process: {:9}\n", nthreads);
     report.reportFullyImplicit(ss);
+
+    if (localsolves_report.total_linearizations > 0) {
+        ss << "======  Accumulated local solve data  ======\n";
+        localsolves_report.reportFullyImplicit(ss);
+    }
+
     OpmLog::info(ss.str());
 }
 

opm/simulators/flow/Banners.hpp

@@ -28,6 +28,7 @@
 namespace Opm {
 
 struct SimulatorReport;
+struct SimulatorReportSingle;
 
 // Print an ASCII-art header to the PRT and DEBUG files.
 void printPRTHeader(const std::string& parameters,
@@ -39,7 +40,8 @@ void printFlowBanner(int nprocs, int threads, std::string_view moduleVersionName
 
 // Print flow application trailer.
 void printFlowTrailer(int nprocs, int nthreads,
-                      const SimulatorReport& report);
+                      const SimulatorReport& report,
+                      const SimulatorReportSingle& localsolves_report);
 
 } // namespace Opm
 

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -177,8 +177,38 @@ template<class TypeTag, class MyTypeTag>
 struct NetworkMaxIterations {
     using type = UndefinedProperty;
 };
-
+template<class TypeTag, class MyTypeTag>
-
+struct NonlinearSolver {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct LocalSolveApproach {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct MaxLocalSolveIterations {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct LocalToleranceScalingMb {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct LocalToleranceScalingCnv {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct NumLocalDomains {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct LocalDomainsPartitioningImbalance {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct LocalDomainsPartitioningMethod {
+    using type = UndefinedProperty;
+};
 template<class TypeTag>
 struct DbhpMaxRel<TypeTag, TTag::FlowModelParameters> {
     using type = GetPropType<TypeTag, Scalar>;
@@ -335,11 +365,42 @@ template<class TypeTag>
 struct NetworkMaxIterations<TypeTag, TTag::FlowModelParameters> {
     static constexpr int value = 200;
 };
-
+template<class TypeTag>
-
+struct NonlinearSolver<TypeTag, TTag::FlowModelParameters> {
-
+    static constexpr auto value = "newton";
-
+};
-
+template<class TypeTag>
+struct LocalSolveApproach<TypeTag, TTag::FlowModelParameters> {
+    static constexpr auto value = "jacobi";
+};
+template<class TypeTag>
+struct MaxLocalSolveIterations<TypeTag, TTag::FlowModelParameters> {
+    static constexpr int value = 20;
+};
+template<class TypeTag>
+struct LocalToleranceScalingMb<TypeTag, TTag::FlowModelParameters> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 1.0;
+};
+template<class TypeTag>
+struct LocalToleranceScalingCnv<TypeTag, TTag::FlowModelParameters> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 0.01;
+};
+template<class TypeTag>
+struct NumLocalDomains<TypeTag, TTag::FlowModelParameters> {
+    using type = int;
+    static constexpr auto value = 0;
+};
+template<class TypeTag>
+struct LocalDomainsPartitioningImbalance<TypeTag, TTag::FlowModelParameters> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr auto value = type{1.03};
+};
+template<class TypeTag>
+struct LocalDomainsPartitioningMethod<TypeTag, TTag::FlowModelParameters> {
+    static constexpr auto value = "zoltan";
+};
 // if openMP is available, determine the number threads per process automatically.
 #if _OPENMP
 template<class TypeTag>
@@ -461,6 +522,19 @@ namespace Opm
         /// Maximum number of iterations in the network solver before giving up
         int network_max_iterations_;
 
+        /// Nonlinear solver type: newton or nldd.
+        std::string nonlinear_solver_;
+        /// 'jacobi' and 'gauss-seidel' supported.
+        std::string local_solve_approach_;
+
+        int max_local_solve_iterations_;
+
+        double local_tolerance_scaling_mb_;
+        double local_tolerance_scaling_cnv_;
+
+        int num_local_domains_{0};
+        double local_domain_partition_imbalance_{1.03};
+        std::string local_domain_partition_method_;
 
         /// Construct from user parameters or defaults.
         BlackoilModelParametersEbos()
@@ -497,6 +571,14 @@ namespace Opm
             check_well_operability_iter_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter);
             max_number_of_well_switches_ = EWOMS_GET_PARAM(TypeTag, int, MaximumNumberOfWellSwitches);
             use_average_density_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseAverageDensityMsWells);
+            nonlinear_solver_ = EWOMS_GET_PARAM(TypeTag, std::string, NonlinearSolver);
+            local_solve_approach_ = EWOMS_GET_PARAM(TypeTag, std::string, LocalSolveApproach);
+            max_local_solve_iterations_ = EWOMS_GET_PARAM(TypeTag, int, MaxLocalSolveIterations);
+            local_tolerance_scaling_mb_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingMb);
+            local_tolerance_scaling_cnv_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingCnv);
+            num_local_domains_ = EWOMS_GET_PARAM(TypeTag, int, NumLocalDomains);
+            local_domain_partition_imbalance_ = std::max(1.0, EWOMS_GET_PARAM(TypeTag, double, LocalDomainsPartitioningImbalance));
+            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);            
@@ -540,6 +622,15 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, bool, UseAverageDensityMsWells, "Approximate segment densitities by averaging over segment and its outlet");
             EWOMS_REGISTER_PARAM(TypeTag, int, NetworkMaxStrictIterations, "Maximum iterations in network solver before relaxing tolerance");
             EWOMS_REGISTER_PARAM(TypeTag, int, NetworkMaxIterations, "Maximum number of iterations in the network solver before giving up");
+            EWOMS_REGISTER_PARAM(TypeTag, std::string, NonlinearSolver, "Choose nonlinear solver. Valid choices are newton or nldd.");
+            EWOMS_REGISTER_PARAM(TypeTag, std::string, LocalSolveApproach, "Choose local solve approach. Valid choices are jacobi and gauss-seidel");
+            EWOMS_REGISTER_PARAM(TypeTag, int, MaxLocalSolveIterations, "Max iterations for local solves with NLDD nonlinear solver.");
+            EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalToleranceScalingMb, "Set lower than 1.0 to use stricter convergence tolerance for local solves.");
+            EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalToleranceScalingCnv, "Set lower than 1.0 to use stricter convergence tolerance for local solves.");
+            EWOMS_REGISTER_PARAM(TypeTag, int, NumLocalDomains, "Number of local domains for NLDD nonlinear solver.");
+            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'.");
         }
     };
 } // namespace Opm

opm/simulators/flow/FlowMainEbos.hpp

@@ -495,7 +495,7 @@ void handleExtraConvergenceOutput(SimulatorReport& report,
                 = omp_get_max_threads();
 #endif
 
-            printFlowTrailer(mpi_size_, threads, report);
+            printFlowTrailer(mpi_size_, threads, report, simulator_->model().localAccumulatedReports());
 
             detail::handleExtraConvergenceOutput(report,
                                                  EWOMS_GET_PARAM(TypeTag, std::string, OutputExtraConvergenceInfo),

opm/simulators/linalg/setupPropertyTree.cpp

@@ -96,6 +96,10 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters
         return setupILU(conf, p);
     }
 
+    if (conf == "umfpack") {
+        return setupUMFPack(conf, p);
+    }
+
     // At this point, the only separate ISAI implementation is with the OpenCL code, and
     // it will check this argument to see if it should be using ISAI. The parameter tree
     // will be ignored, so this is just a dummy configuration to avoid the throw below.
@@ -263,4 +267,15 @@ setupILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParamet
 }
 
 
+PropertyTree
+setupUMFPack([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
+{
+    using namespace std::string_literals;
+    PropertyTree prm;
+    prm.put("verbosity", p.linear_solver_verbosity_);
+    prm.put("solver", "umfpack"s);
+    return prm;
+}
+
+
 } // namespace Opm

opm/simulators/linalg/setupPropertyTree.hpp

@@ -37,6 +37,7 @@ PropertyTree setupCPRW(const std::string& conf, const FlowLinearSolverParameters
 PropertyTree setupCPR(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupAMG(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupILU(const std::string& conf, const FlowLinearSolverParameters& p);
+PropertyTree setupUMFPack(const std::string& conf, const FlowLinearSolverParameters& p);
 
 } // namespace Opm
 

opm/simulators/timestepping/SimulatorReport.cpp

@@ -240,6 +240,7 @@ namespace Opm
 
     void SimulatorReport::reportFullyImplicit(std::ostream& os) const
     {
+        os << fmt::format("Number of timesteps:     {:9}\n", stepreports.size());
         success.reportFullyImplicit(os, &failure);
     }
 

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1049,7 +1049,7 @@ namespace Opm {
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
-    assembleDomain(const int /*iterationIdx*/,
+    assembleDomain([[maybe_unused]] const int iterationIdx,
                    const double dt,
                    const Domain& domain)
     {
@@ -1575,7 +1575,6 @@ namespace Opm {
         const bool relax_tolerance = iterationIdx > param_.strict_outer_iter_wells_;
 
         Opm::DeferredLogger local_deferredLogger;
-
         ConvergenceReport local_report;
         for (const auto& well : well_container_) {
             if ((well_domain_.at(well->name()) == domain.index)) {