From a2b96eaa8fc26a144c413749f9b41762e441482e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 9 Jun 2023 13:46:45 +0200 Subject: [PATCH] Clean up linear solver parameters and ISTLSolverEbos. --- .../linalg/FlowLinearSolverParameters.hpp | 99 ++++++++++--------- opm/simulators/linalg/ISTLSolverEbos.hpp | 53 +++++++--- 2 files changed, 92 insertions(+), 60 deletions(-) diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.hpp b/opm/simulators/linalg/FlowLinearSolverParameters.hpp index 425c41a2c..055c2882c 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.hpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.hpp @@ -51,10 +51,6 @@ struct RelaxedLinearSolverReduction { using type = UndefinedProperty; }; -template -struct IluRelaxation { - using type = UndefinedProperty; -}; template struct LinearSolverMaxIter { using type = UndefinedProperty; @@ -67,6 +63,10 @@ struct LinearSolverRestart { // LinearSolverVerbosity defined in opm-models // template +struct IluRelaxation { + using type = UndefinedProperty; +}; +template struct IluFillinLevel { using type = UndefinedProperty; }; @@ -91,11 +91,15 @@ struct LinearSolverIgnoreConvergenceFailure{ using type = UndefinedProperty; }; template -struct PreconditionerAddWellContributions { +struct ScaleLinearSystem { using type = UndefinedProperty; }; template -struct ScaleLinearSystem { +struct LinearSolver { + using type = UndefinedProperty; +}; +template +struct LinearSolverPrintJsonDefinition { using type = UndefinedProperty; }; template @@ -107,10 +111,6 @@ struct CprReuseInterval { using type = UndefinedProperty; }; template -struct LinearSolver { - using type = UndefinedProperty; -}; -template struct AcceleratorMode { using type = UndefinedProperty; }; @@ -137,11 +137,6 @@ struct RelaxedLinearSolverReduction { static constexpr type value = 1e-2; }; template -struct IluRelaxation { - using type = GetPropType; - static constexpr type value = 0.9; -}; -template struct LinearSolverMaxIter { static constexpr int value = 200; }; @@ -154,6 +149,11 @@ struct LinearSolverVerbosity { static constexpr int value = 0; }; template +struct IluRelaxation { + using type = GetPropType; + static constexpr type value = 0.9; +}; +template struct IluFillinLevel { static constexpr int value = 0; }; @@ -178,18 +178,18 @@ struct LinearSolverIgnoreConvergenceFailure static constexpr bool value = false; }; template -struct LinearSolverBackend { - using type = ISTLSolverEbos; -}; -template -struct PreconditionerAddWellContributions { - static constexpr bool value = false; -}; -template struct ScaleLinearSystem { static constexpr bool value = false; }; template +struct LinearSolver { + static constexpr auto value = "ilu0"; +}; +template +struct LinearSolverPrintJsonDefinition { + static constexpr auto value = true; +}; +template struct CprReuseSetup { static constexpr int value = 4; }; @@ -198,10 +198,6 @@ struct CprReuseInterval { static constexpr int value = 30; }; template -struct LinearSolver { - static constexpr auto value = "ilu0"; -}; -template struct AcceleratorMode { static constexpr auto value = "none"; }; @@ -218,6 +214,12 @@ struct OpenclIluParallel { static constexpr bool value = true; // note: false should only be used in debug }; +// Set the backend to be used. +template +struct LinearSolverBackend { + using type = ISTLSolverEbos; +}; + } // namespace Opm::Properties namespace Opm @@ -229,24 +231,24 @@ namespace Opm { double linear_solver_reduction_; double relaxed_linear_solver_reduction_; - double ilu_relaxation_; int linear_solver_maxiter_; int linear_solver_restart_; int linear_solver_verbosity_; + double ilu_relaxation_; int ilu_fillin_level_; MILU_VARIANT ilu_milu_; bool ilu_redblack_; bool ilu_reorder_sphere_; bool newton_use_gmres_; - bool require_full_sparsity_pattern_; bool ignoreConvergenceFailure_; bool scale_linear_system_; std::string linsolver_; + bool linear_solver_print_json_definition_; + int cpr_reuse_setup_; + int cpr_reuse_interval_; std::string accelerator_mode_; int bda_device_id_; int opencl_platform_id_; - int cpr_reuse_setup_; - int cpr_reuse_interval_; bool opencl_ilu_parallel_; template @@ -255,10 +257,10 @@ namespace Opm // TODO: these parameters have undocumented non-trivial dependencies linear_solver_reduction_ = EWOMS_GET_PARAM(TypeTag, double, LinearSolverReduction); relaxed_linear_solver_reduction_ = EWOMS_GET_PARAM(TypeTag, double, RelaxedLinearSolverReduction); - ilu_relaxation_ = EWOMS_GET_PARAM(TypeTag, double, IluRelaxation); linear_solver_maxiter_ = EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIter); linear_solver_restart_ = EWOMS_GET_PARAM(TypeTag, int, LinearSolverRestart); linear_solver_verbosity_ = EWOMS_GET_PARAM(TypeTag, int, LinearSolverVerbosity); + ilu_relaxation_ = EWOMS_GET_PARAM(TypeTag, double, IluRelaxation); ilu_fillin_level_ = EWOMS_GET_PARAM(TypeTag, int, IluFillinLevel); ilu_milu_ = convertString2Milu(EWOMS_GET_PARAM(TypeTag, std::string, MiluVariant)); ilu_redblack_ = EWOMS_GET_PARAM(TypeTag, bool, IluRedblack); @@ -266,10 +268,12 @@ namespace Opm newton_use_gmres_ = EWOMS_GET_PARAM(TypeTag, bool, UseGmres); ignoreConvergenceFailure_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure); scale_linear_system_ = EWOMS_GET_PARAM(TypeTag, bool, ScaleLinearSystem); + linsolver_ = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolver); + linear_solver_print_json_definition_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverPrintJsonDefinition); cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup); cpr_reuse_interval_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseInterval); - if (!EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter) && cprRequestedInDataFile) { + if (!EWOMS_PARAM_IS_SET(TypeTag, std::string, LinearSolver) && cprRequestedInDataFile) { linsolver_ = "cpr"; } else { linsolver_ = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolver); @@ -284,12 +288,12 @@ namespace Opm template static void registerParameters() { - EWOMS_REGISTER_PARAM(TypeTag, double, LinearSolverReduction, "The minimum reduction of the residual which the linear solver must achieve should try to achive"); + EWOMS_REGISTER_PARAM(TypeTag, double, LinearSolverReduction, "The minimum reduction of the residual which the linear solver must achieve"); EWOMS_REGISTER_PARAM(TypeTag, double, RelaxedLinearSolverReduction, "The minimum reduction of the residual which the linear solver need to achieve for the solution to be accepted"); - EWOMS_REGISTER_PARAM(TypeTag, double, IluRelaxation, "The relaxation factor of the linear solver's ILU preconditioner"); EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverMaxIter, "The maximum number of iterations of the linear solver"); EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverRestart, "The number of iterations after which GMRES is restarted"); EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverVerbosity, "The verbosity level of the linear solver (0: off, 2: all)"); + EWOMS_REGISTER_PARAM(TypeTag, double, IluRelaxation, "The relaxation factor of the linear solver's ILU preconditioner"); EWOMS_REGISTER_PARAM(TypeTag, int, IluFillinLevel, "The fill-in level of the linear solver's ILU preconditioner"); EWOMS_REGISTER_PARAM(TypeTag, std::string, MiluVariant, "Specify which variant of the modified-ILU preconditioner ought to be used. Possible variants are: ILU (default, plain ILU), MILU_1 (lump diagonal with dropped row entries), MILU_2 (lump diagonal with the sum of the absolute values of the dropped row entries), MILU_3 (if diagonal is positive add sum of dropped row entrires. Otherwise subtract them), MILU_4 (if diagonal is positive add sum of dropped row entrires. Otherwise do nothing"); EWOMS_REGISTER_PARAM(TypeTag, bool, IluRedblack, "Use red-black partitioning for the ILU preconditioner"); @@ -297,9 +301,10 @@ 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, 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."); - 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, AcceleratorMode, "Choose a linear solver, usage: '--accelerator-mode=[none|cusparse|opencl|amgcl|rocalution]'"); EWOMS_REGISTER_PARAM(TypeTag, int, BdaDeviceId, "Choose device ID for cusparseSolver or openclSolver, use 'nvidia-smi' or 'clinfo' to determine valid IDs"); EWOMS_REGISTER_PARAM(TypeTag, int, OpenclPlatformId, "Choose platform ID for openclSolver, use 'clinfo' to determine valid platform IDs"); @@ -311,19 +316,23 @@ namespace Opm // set default values void reset() { - newton_use_gmres_ = false; - linear_solver_reduction_ = 1e-2; relaxed_linear_solver_reduction_ = 1e-2; - linear_solver_maxiter_ = 150; - linear_solver_restart_ = 40; - linear_solver_verbosity_ = 0; - require_full_sparsity_pattern_ = false; - ignoreConvergenceFailure_ = false; - ilu_fillin_level_ = 0; + linear_solver_reduction_ = 1e-2; + linear_solver_maxiter_ = 200; + linear_solver_restart_ = 40; + linear_solver_verbosity_ = 0; ilu_relaxation_ = 0.9; + ilu_fillin_level_ = 0; ilu_milu_ = MILU_VARIANT::ILU; ilu_redblack_ = false; - ilu_reorder_sphere_ = true; + ilu_reorder_sphere_ = false; + newton_use_gmres_ = false; + ignoreConvergenceFailure_ = false; + scale_linear_system_ = false; + linsolver_ = "ilu0"; + linear_solver_print_json_definition_ = true; + cpr_reuse_setup_ = 4; + cpr_reuse_interval_ = 30; accelerator_mode_ = "none"; bda_device_id_ = 0; opencl_platform_id_ = 0; diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 068a3ce93..f4d9194ea 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -231,25 +231,45 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, } /// Construct a system solver. - /// \param[in] parallelInformation In the case of a parallel run - /// with dune-istl the information about the parallelization. + /// \param[in] simulator The opm-models simulator object + /// \param[in] parameters Explicit parameters for solver setup, do not + /// read them from command line parameters. + ISTLSolverEbos(const Simulator& simulator, const FlowLinearSolverParameters& parameters) + : simulator_(simulator), + iterations_( 0 ), + calls_( 0 ), + converged_(false), + matrix_(nullptr), + parameters_(parameters) + { + initialize(); + } + + /// Construct a system solver. + /// \param[in] simulator The opm-models simulator object explicit ISTLSolverEbos(const Simulator& simulator) : simulator_(simulator), iterations_( 0 ), calls_( 0 ), converged_(false), - matrix_() + matrix_(nullptr) + { + parameters_.template init(simulator_.vanguard().eclState().getSimulationConfig().useCPR()); + initialize(); + } + + void initialize() { OPM_TIMEBLOCK(IstlSolverEbos); - const bool on_io_rank = (simulator.gridView().comm().rank() == 0); -#if HAVE_MPI - comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) ); -#endif - parameters_.template init(simulator_.vanguard().eclState().getSimulationConfig().useCPR()); prm_ = setupPropertyTree(parameters_, EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter), EWOMS_PARAM_IS_SET(TypeTag, double, LinearSolverReduction)); + const bool on_io_rank = (simulator_.gridView().comm().rank() == 0); +#if HAVE_MPI + comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) ); +#endif + #if COMPILE_BDA_BRIDGE { std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode); @@ -300,7 +320,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, flexibleSolver_.interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true); // Print parameters to PRT/DBG logs. - if (on_io_rank) { + if (on_io_rank && parameters_.linear_solver_print_json_definition_) { std::ostringstream os; os << "Property tree for linear solver:\n"; prm_.write_json(os, true); @@ -314,12 +334,17 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, } void prepare(const SparseMatrixAdapter& M, Vector& b) + { + prepare(M.istlMatrix(), b); + } + + void prepare(const Matrix& M, Vector& b) { OPM_TIMEBLOCK(istlSolverEbosPrepare); - static bool firstcall = true; + const bool firstcall = (matrix_ == nullptr); #if HAVE_MPI if (firstcall) { - const size_t size = M.istlMatrix().N(); + const size_t size = M.N(); detail::copyParValues(parallelInformation_, size, *comm_); } #endif @@ -329,7 +354,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, // ebos will not change the matrix object. Hence simply store a pointer // to the original one with a deleter that does nothing. // Outch! We need to be able to scale the linear system! Hence const_cast - matrix_ = const_cast(&M.istlMatrix()); + matrix_ = const_cast(&M); useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver) @@ -343,7 +368,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, #endif } else { // Pointers should not change - if ( &(M.istlMatrix()) != matrix_ ) { + if ( &M != matrix_ ) { OPM_THROW(std::logic_error, "Matrix objects are expected to be reused when reassembling!"); } @@ -360,8 +385,6 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, #else prepareFlexibleSolver(); #endif - - firstcall = false; }