From cede2b4add66e56ab23a2bcd45bbdd39a60aa083 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 31 Aug 2023 10:04:18 +0200 Subject: [PATCH] Improving generality and output of hybrid approach. --- opm/simulators/flow/BlackoilModelEbos.hpp | 76 +++++++------------ .../linalg/FlowLinearSolverParameters.hpp | 2 +- opm/simulators/linalg/ISTLSolverEbos.hpp | 15 ++-- 3 files changed, 35 insertions(+), 58 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index cebe916f1..de556497f 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/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 times(numSolvers); + std::vector setupTimes(numSolvers); - auto ebosJacCopy = ebosJac; - auto ebosResidCopy = ebosResid; - - // solve with linear solver 0 x = 0.0; - BVector x0(x); - ebosSolver.setActiveSolver(0); - perfTimer.start(); - ebosSolver.prepare(ebosJac, ebosResid); - auto linear_solve_setup_time0 = perfTimer.stop(); - perfTimer.reset(); - ebosSolver.setResidual(ebosResid); - perfTimer.start(); - ebosSolver.solve(x0); - auto time0 = perfTimer.stop(); - perfTimer.reset(); - - // reset matrix to the original and Rhs - ebosJac = ebosJacCopy; - ebosResid = ebosResidCopy; - - // 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"); + std::vector x_trial(numSolvers, x); + for (int solver = 0; solver < numSolvers; ++solver) { + BVector x0(x); + ebosSolver.setActiveSolver(solver); + perfTimer.start(); + ebosSolver.prepare(ebosJac, ebosResid); + setupTimes[solver] = perfTimer.stop(); + perfTimer.reset(); + ebosSolver.setResidual(ebosResid); + perfTimer.start(); + ebosSolver.solve(x_trial[solver]); + times[solver] = perfTimer.stop(); + perfTimer.reset(); + if (terminal_output_) { + OpmLog::debug(fmt::format("Solver time {}: {}", solver, times[solver])); + } } - 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_ ) { - OpmLog::debug("Using solver " + std::to_string(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; - } + linear_solve_setup_time_ = setupTimes[fastest_solver]; + x = x_trial[fastest_solver]; ebosSolver.setActiveSolver(fastest_solver); } else { diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.hpp b/opm/simulators/linalg/FlowLinearSolverParameters.hpp index 651f3c1c4..eeacbfd4d 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.hpp +++ b/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."); diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 6628930a5..dc20031f2 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -211,12 +211,12 @@ std::unique_ptr 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 blockJacobiAdjacency(const Grid& grid, void setActiveSolver(const int num) { - if (num>prm_.size()-1) { - OPM_THROW(std::logic_error,"Solver number"+std::to_string(num)+"not available."); + if (num > static_cast(prm_.size()) - 1) { + 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_ + ")"); } }