Improving generality and output of hybrid approach.

opm/simulators/flow/BlackoilModelEbos.hpp

@@ -596,65 +596,41 @@ namespace Opm {
             auto& ebosResid = ebosSimulator_.model().linearizer().residual();
             auto& ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver();
 
-            if ((ebosSolver.numAvailableSolvers() > 1) && (ebosSolver.getSolveCount() % 100 == 0)) {
+            const int numSolvers = ebosSolver.numAvailableSolvers();
+            if ((numSolvers > 1) && (ebosSolver.getSolveCount() % 100 == 0)) {
 
                 if ( terminal_output_ ) {
-                    OpmLog::debug("Running speed test for comparing available linear solvers.");
+                    OpmLog::debug("\nRunning speed test for comparing available linear solvers.");
                 }
+
                 Dune::Timer perfTimer;
+                std::vector<double> times(numSolvers);
+                std::vector<double> setupTimes(numSolvers);
 
-                auto ebosJacCopy = ebosJac;
-                auto ebosResidCopy = ebosResid;
-
-                // solve with linear solver 0
                 x = 0.0;
-                BVector x0(x);
+                std::vector<BVector> x_trial(numSolvers, x);
-                ebosSolver.setActiveSolver(0);
+                for (int solver = 0; solver < numSolvers; ++solver) {
-                perfTimer.start();
+                    BVector x0(x);
-                ebosSolver.prepare(ebosJac, ebosResid);
+                    ebosSolver.setActiveSolver(solver);
-                auto linear_solve_setup_time0 = perfTimer.stop();
+                    perfTimer.start();
-                perfTimer.reset();
+                    ebosSolver.prepare(ebosJac, ebosResid);
-                ebosSolver.setResidual(ebosResid);
+                    setupTimes[solver] = perfTimer.stop();
-                perfTimer.start();
+                    perfTimer.reset();
-                ebosSolver.solve(x0);
+                    ebosSolver.setResidual(ebosResid);
-                auto time0 = perfTimer.stop();
+                    perfTimer.start();
-                perfTimer.reset();
+                    ebosSolver.solve(x_trial[solver]);
-
+                    times[solver] = perfTimer.stop();
-                // reset matrix to the original and Rhs
+                    perfTimer.reset();
-                ebosJac = ebosJacCopy;
+                    if (terminal_output_) {
-                ebosResid = ebosResidCopy;
+                        OpmLog::debug(fmt::format("Solver time {}: {}", solver, times[solver]));
-
+                    }
-                // solve with linear solver 1
-                x = 0.0;
-                BVector x1(x);
-                ebosSolver.setActiveSolver(1);
-                perfTimer.start();
-                ebosSolver.prepare(ebosJac, ebosResid);
-                auto linear_solve_setup_time1 = perfTimer.stop();
-                perfTimer.reset();
-                ebosSolver.setResidual(ebosResid);
-                perfTimer.start();
-                ebosSolver.solve(x1);
-                auto time1 = perfTimer.stop();
-                perfTimer.reset();
-
-                if ( terminal_output_ ) {
-                    OpmLog::debug("Solver time 0: " + std::to_string(time0) + "\n");
-                    OpmLog::debug("Solver time 1: " + std::to_string(time1) + "\n");
                 }
 
-                int fastest_solver = time0 < time1 ? 0 : 1;
+                int fastest_solver = std::min_element(times.begin(), times.end()) - times.begin();
+                // Use timing on rank 0 to determine fastest, must be consistent across ranks.
                 grid_.comm().broadcast(&fastest_solver, 1, 0);
-                if ( terminal_output_ ) {
+                linear_solve_setup_time_ = setupTimes[fastest_solver];
-                    OpmLog::debug("Using solver " + std::to_string(fastest_solver));
+                x = x_trial[fastest_solver];
-                }
-                if (fastest_solver == 0) {
-                    linear_solve_setup_time_ = linear_solve_setup_time0;
-                    x=x0;
-                } else {
-                    linear_solve_setup_time_ = linear_solve_setup_time1;
-                    x=x1;
-                }
                 ebosSolver.setActiveSolver(fastest_solver);
 
             } else {

opm/simulators/linalg/FlowLinearSolverParameters.hpp

@@ -307,7 +307,7 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, bool, UseGmres, "Use GMRES as the linear solver");
             EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge");
             EWOMS_REGISTER_PARAM(TypeTag, bool, ScaleLinearSystem, "Scale linear system according to equation scale and primary variable types");
-            EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolver, "Configuration of solver. Valid options are: ilu0 (default), cprw, cpr (an alias for cprw), cpr_quasiimpes, cpr_trueimpes or amg. Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'");
+            EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolver, "Configuration of solver. Valid options are: ilu0 (default), cprw, cpr (an alias for cprw), cpr_quasiimpes, cpr_trueimpes, amg or hybrid (experimental). Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'");
             EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverPrintJsonDefinition, "Write the JSON definition of the linear solver setup to the DBG file.");
             EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse preconditioner setup. Valid options are 0: recreate the preconditioner for every linear solve, 1: recreate once every timestep, 2: recreate if last linear solve took more than 10 iterations, 3: never recreate, 4: recreated every CprReuseInterval");
             EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseInterval, "Reuse preconditioner interval. Used when CprReuseSetup is set to 4, then the preconditioner will be fully recreated instead of reused every N linear solve, where N is this parameter.");

opm/simulators/linalg/ISTLSolverEbos.hpp

@@ -211,12 +211,12 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
         void initialize(bool have_gpu = false)
         {
             OPM_TIMEBLOCK(IstlSolverEbos);
-            // prm_ = setupPropertyTree(parameters_,
-            //                          EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter),
-            //                          EWOMS_PARAM_IS_SET(TypeTag, double, LinearSolverReduction));
 
             if (parameters_[0].linsolver_ == "hybrid") {
-                // ------------ the following is hard coded for testing purposes!!!
+                // Experimental hybrid configuration.
+                // When chosen, will set up two solvers, one with CPRW
+                // and the other with ILU0 preconditioner. More general
+                // options may be added later.
                 prm_.clear();
                 parameters_.clear();
                 {
@@ -299,13 +299,14 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
 
         void setActiveSolver(const int num)
         {
-            if (num>prm_.size()-1) {
+            if (num > static_cast<int>(prm_.size()) - 1) {
-                OPM_THROW(std::logic_error,"Solver number"+std::to_string(num)+"not available.");
+                OPM_THROW(std::logic_error, "Solver number " + std::to_string(num) + " not available.");
             }
             activeSolverNum_ = num;
             auto cc = Dune::MPIHelper::getCollectiveCommunication();
             if (cc.rank() == 0) {
-                OpmLog::note("Active solver = "+std::to_string(activeSolverNum_)+"\n");
+                OpmLog::debug("Active solver = " + std::to_string(activeSolverNum_)
+                              + " (" + parameters_[activeSolverNum_].linsolver_ + ")");
             }
         }