OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4700 from atgeirr/linear-solver-improvements

Browse Source 
Linear solver improvements
This commit is contained in:
Atgeirr Flø Rasmussen 2023-06-11 12:43:36 +02:00 committed by GitHub
commit d926f17abe
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 122 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
opm/simulators/linalg/FlexibleSolver.hpp
View File

@ -88,6 +88,8 @@ private:
void initSolver(const Opm::PropertyTree& prm, const bool is_iorank);
void recreateDirectSolver();
// Main initialization routine.
// Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
template <class Comm>
@ -101,6 +103,7 @@ private:
std::shared_ptr<AbstractPrecondType> preconditioner_;
std::shared_ptr<AbstractScalarProductType> scalarproduct_;
std::shared_ptr<AbstractSolverType> linsolver_;
bool direct_solver_ = false;
};
} // namespace Dune

29
opm/simulators/linalg/FlexibleSolver_impl.hpp
View File

@ -73,6 +73,9 @@ namespace Dune
FlexibleSolver<Operator>::
apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res)
{
if (direct_solver_) {
recreateDirectSolver();
}
linsolver_->apply(x, rhs, res);
}
@ -81,6 +84,9 @@ namespace Dune
FlexibleSolver<Operator>::
apply(VectorType& x, VectorType& rhs, double reduction, Dune::InverseOperatorResult& res)
{
if (direct_solver_) {
recreateDirectSolver();
}
linsolver_->apply(x, rhs, reduction, res);
}
@ -185,9 +191,9 @@ namespace Dune
verbosity);
#if HAVE_SUITESPARSE_UMFPACK
} else if (solver_type == "umfpack") {
bool dummy = false;
using MatrixType = std::remove_const_t<std::remove_reference_t<decltype(linearoperator_for_solver_->getmat())>>;
linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), verbosity, dummy);
linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), verbosity, false);
direct_solver_ = true;
#endif
#if HAVE_CUDA
} else if (solver_type == "cubicgstab") {
@ -206,6 +212,25 @@ namespace Dune
}
// For now, the only direct solver we support is UMFPACK from SuiteSparse.
// When the matrix is updated (keeping sparsity pattern) it is possible to
// exploit separation of symbolic and numeric factorization, but we do not
// do so at this point. For complete generality, the solver abstract class
// Dune::InverseOperator<> should be extended with an update() function.
template <class Operator>
void
FlexibleSolver<Operator>::
recreateDirectSolver()
{
#if HAVE_SUITESPARSE_UMFPACK
using MatrixType = std::remove_const_t<std::remove_reference_t<decltype(linearoperator_for_solver_->getmat())>>;
linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), 0, false);
#else
OPM_THROW(std::logic_error, "Direct solver specified, but the FlexibleSolver class was not compiled with SuiteSparse support.");
#endif
}
// Main initialization routine.
// Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
template <class Operator>

99
opm/simulators/linalg/FlowLinearSolverParameters.hpp
View File

@ -51,10 +51,6 @@ struct RelaxedLinearSolverReduction {
using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct IluRelaxation {
using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct LinearSolverMaxIter {
using type = UndefinedProperty;
@ -67,6 +63,10 @@ struct LinearSolverRestart {
// LinearSolverVerbosity defined in opm-models
//
template<class TypeTag, class MyTypeTag>
struct IluRelaxation {
using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct IluFillinLevel {
using type = UndefinedProperty;
};
@ -91,11 +91,15 @@ struct LinearSolverIgnoreConvergenceFailure{
using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct PreconditionerAddWellContributions {
struct ScaleLinearSystem {
using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct ScaleLinearSystem {
struct LinearSolver {
using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct LinearSolverPrintJsonDefinition {
using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
@ -107,10 +111,6 @@ struct CprReuseInterval {
using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct LinearSolver {
using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct AcceleratorMode {
using type = UndefinedProperty;
};
@ -137,11 +137,6 @@ struct RelaxedLinearSolverReduction<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr type value = 1e-2;
};
template<class TypeTag>
struct IluRelaxation<TypeTag, TTag::FlowIstlSolverParams> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 0.9;
};
template<class TypeTag>
struct LinearSolverMaxIter<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr int value = 200;
};
@ -154,6 +149,11 @@ struct LinearSolverVerbosity<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr int value = 0;
};
template<class TypeTag>
struct IluRelaxation<TypeTag, TTag::FlowIstlSolverParams> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 0.9;
};
template<class TypeTag>
struct IluFillinLevel<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr int value = 0;
};
@ -178,18 +178,18 @@ struct LinearSolverIgnoreConvergenceFailure<TypeTag, TTag::FlowIstlSolverParams>
static constexpr bool value = false;
};
template<class TypeTag>
struct LinearSolverBackend<TypeTag, TTag::FlowIstlSolverParams> {
using type = ISTLSolverEbos<TypeTag>;
};
template<class TypeTag>
struct PreconditionerAddWellContributions<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr bool value = false;
};
template<class TypeTag>
struct ScaleLinearSystem<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr bool value = false;
};
template<class TypeTag>
struct LinearSolver<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr auto value = "ilu0";
};
template<class TypeTag>
struct LinearSolverPrintJsonDefinition<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr auto value = true;
};
template<class TypeTag>
struct CprReuseSetup<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr int value = 4;
};
@ -198,10 +198,6 @@ struct CprReuseInterval<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr int value = 30;
};
template<class TypeTag>
struct LinearSolver<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr auto value = "ilu0";
};
template<class TypeTag>
struct AcceleratorMode<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr auto value = "none";
};
@ -218,6 +214,12 @@ struct OpenclIluParallel<TypeTag, TTag::FlowIstlSolverParams> {
static constexpr bool value = true; // note: false should only be used in debug
};
// Set the backend to be used.
template<class TypeTag>
struct LinearSolverBackend<TypeTag, TTag::FlowIstlSolverParams> {
using type = ISTLSolverEbos<TypeTag>;
};
} // 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 <class TypeTag>
@ -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 <class TypeTag>
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;

53
opm/simulators/linalg/ISTLSolverEbos.hpp
View File

@ -231,25 +231,45 @@ std::unique_ptr<Matrix> 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<TypeTag>(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<TypeTag>(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<Matrix> 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<Matrix> 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<Matrix> 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<Matrix*>(&M.istlMatrix());
matrix_ = const_cast<Matrix*>(&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<Matrix> 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<Matrix> blockJacobiAdjacency(const Grid& grid,
#else
prepareFlexibleSolver();
#endif
firstcall = false;
}