diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index e643f9ad4..ee18a087c 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -77,9 +77,8 @@ Dune::MatrixAdapter createOperator(const Dune::MatrixAdapter&, con template std::unique_ptr< Dune::MatrixAdapter > createOperatorPtr(const Dune::MatrixAdapter&, const M& matrix, const T&) { - return std::make_unique< Dune::MatrixAdapter >(matrix); -} - + return std::make_unique< Dune::MatrixAdapter >(matrix); +} /** * \brief Creates an OverlappingSchwarzOperator as an operator. * @@ -119,13 +118,9 @@ scaleMatrixDRS(const Operator& op, const Communication& comm, for (auto j = (*i).begin(); j != endj; ++j) { BlockVector bvec(0.0); Block& block = *j; - for (std::size_t ii = 0; ii < Block::rows; ii++) { - for (std::size_t jj = 0; jj < Block::cols; jj++) { - // should introduce limmits which also change the weights - bvec[jj] += bw[ii]*block[ii][jj]; - } - } - block[pressureIndex] = bvec; + BlockVector& bvec = block[pressureIndex]; + // should introduce limits which also change the weights + block.mtv(bw, bvec); } } } @@ -144,13 +139,9 @@ void scaleVectorDRS(Vector& vector, std::size_t pressureIndex, const Opm::CPRPar using Block = typename Vector::block_type; if (param.cpr_use_drs_) { for (std::size_t j = 0; j < vector.size(); ++j) { - double val(0.0); Block& block = vector[j]; const Block& bw = weights[j]; - for (std::size_t i = 0; i < Block::dimension; i++) { - val += bw[i]*block[i]; - } - block[pressureIndex] = val; + block[pressureIndex] = bw.dot(block); } } } @@ -468,8 +459,8 @@ private: } #endif } - else - { + else + { #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) Dune::LoopSolver solver(const_cast(op_), *sp, *prec, tolerance, maxit, verbosity); @@ -489,9 +480,9 @@ private: tolerance, maxit, verbosity); solver.apply(x,b,res); } - -#endif - } + +#endif + } #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) delete sp; @@ -983,20 +974,19 @@ public: * \param comm The information about the parallelization. */ BlackoilAmg(const CPRParameter& param, - const typename TwoLevelMethod::FineDomainType& weights, + const typename TwoLevelMethod::FineDomainType& weights, const Operator& fineOperator, const Criterion& criterion, const SmootherArgs& smargs, const Communication& comm) : param_(param), - weights_(weights), + weights_(weights), scaledMatrixOperator_(Detail::scaleMatrixDRS(fineOperator, comm, - COMPONENT_INDEX, weights, param)), + COMPONENT_INDEX, weights, param)), smoother_(Detail::constructSmoother(std::get<1>(scaledMatrixOperator_), - smargs, comm)), + smargs, comm)), levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_), coarseSolverPolicy_(¶m, smargs, criterion), twoLevelMethod_(std::get<1>(scaledMatrixOperator_), smoother_, - levelTransferPolicy_, - coarseSolverPolicy_, 0, 1) + levelTransferPolicy_, coarseSolverPolicy_, 0, 1) {} void pre(typename TwoLevelMethod::FineDomainType& x, diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index 81dbff274..1e46558f0 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -188,7 +188,7 @@ inline void createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::unique_ptr< BlackoilAmg >& amgPtr, const CPRParameter& params, - const Vector& weights) + const Vector& weights) { using AMG = BlackoilAmg; const int verbosity = ( params.cpr_solver_verbose_ && comm.communicator().rank() == 0 ) ? 1 : 0; diff --git a/opm/autodiff/FlowLinearSolverParameters.hpp b/opm/autodiff/FlowLinearSolverParameters.hpp index 7ab940cf8..feeb77d7f 100644 --- a/opm/autodiff/FlowLinearSolverParameters.hpp +++ b/opm/autodiff/FlowLinearSolverParameters.hpp @@ -128,14 +128,14 @@ namespace Opm cpr_ilu_redblack_ = false; cpr_ilu_reorder_sphere_ = true; cpr_max_ell_iter_ = 25; - cpr_ell_solvetype_ = 0; - cpr_use_drs_ = false; - cpr_max_iter_ = 25; + cpr_ell_solvetype_ = 0; + cpr_use_drs_ = false; + cpr_max_iter_ = 25; cpr_use_amg_ = true; cpr_use_bicgstab_ = true; cpr_solver_verbose_ = 0; cpr_pressure_aggregation_ = false; - cpr_reuse_setup_ = 0; + cpr_reuse_setup_ = 0; } }; @@ -181,14 +181,13 @@ namespace Opm ignoreConvergenceFailure_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure); linear_solver_use_amg_ = EWOMS_GET_PARAM(TypeTag, bool, UseAmg); use_cpr_ = EWOMS_GET_PARAM(TypeTag, bool, UseCpr); - system_strategy_ = EWOMS_GET_PARAM(TypeTag, std::string, SystemStrategy); - scale_linear_system_ = EWOMS_GET_PARAM(TypeTag, bool, ScaleLinearSystem); - cpr_solver_verbose_ = EWOMS_GET_PARAM(TypeTag, int, CprSolverVerbose); - cpr_use_drs_ = EWOMS_GET_PARAM(TypeTag, bool, CprUseDrs); - cpr_max_iter_ = EWOMS_GET_PARAM(TypeTag, int, CprMaxIter); - cpr_ell_solvetype_ = EWOMS_GET_PARAM(TypeTag, int, CprEllSolvetype); - cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup); - + system_strategy_ = EWOMS_GET_PARAM(TypeTag, std::string, SystemStrategy); + scale_linear_system_ = EWOMS_GET_PARAM(TypeTag, bool, ScaleLinearSystem); + cpr_solver_verbose_ = EWOMS_GET_PARAM(TypeTag, int, CprSolverVerbose); + cpr_use_drs_ = EWOMS_GET_PARAM(TypeTag, bool, CprUseDrs); + cpr_max_iter_ = EWOMS_GET_PARAM(TypeTag, int, CprMaxIter); + cpr_ell_solvetype_ = EWOMS_GET_PARAM(TypeTag, int, CprEllSolvetype); + cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup); } template @@ -208,14 +207,13 @@ namespace Opm 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, UseAmg, "Use AMG as the linear solver's preconditioner"); EWOMS_REGISTER_PARAM(TypeTag, bool, UseCpr, "Use CPR as the linear solver's preconditioner"); - EWOMS_REGISTER_PARAM(TypeTag, std::string, SystemStrategy, "Strategy for reformulating and scale linear system"); - EWOMS_REGISTER_PARAM(TypeTag, bool, ScaleLinearSystem, "Scale linear system according to equation scale and primary variable types"); - EWOMS_REGISTER_PARAM(TypeTag, int, CprSolverVerbose, "Verbose for cpr solver"); - EWOMS_REGISTER_PARAM(TypeTag, bool, CprUseDrs, "Use dynamic row sum using weighs"); - EWOMS_REGISTER_PARAM(TypeTag, int, CprMaxIter, "MaxIterations of the pressure amg solver"); - EWOMS_REGISTER_PARAM(TypeTag, int, CprEllSolvetype, "solver type of elliptic solve 0 bicgstab 1 cg other only amg preconditioner"); - EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse Amg Setup"); - + EWOMS_REGISTER_PARAM(TypeTag, std::string, SystemStrategy, "Strategy for reformulating and scale linear system"); + EWOMS_REGISTER_PARAM(TypeTag, bool, ScaleLinearSystem, "Scale linear system according to equation scale and primary variable types"); + EWOMS_REGISTER_PARAM(TypeTag, int, CprSolverVerbose, "Verbose for cpr solver"); + EWOMS_REGISTER_PARAM(TypeTag, bool, CprUseDrs, "Use dynamic row sum using weighs"); + EWOMS_REGISTER_PARAM(TypeTag, int, CprMaxIter, "MaxIterations of the pressure amg solver"); + EWOMS_REGISTER_PARAM(TypeTag, int, CprEllSolvetype, "solver type of elliptic solve 0 bicgstab 1 cg other only amg preconditioner"); + EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse Amg Setup"); } FlowLinearSolverParameters() { reset(); } diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index 9c10b7f50..d59a02b88 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -253,6 +253,10 @@ protected: parameters_.cpr_use_drs_ = false; } } else { + if (parameters_.use_cpr_ && parameters_.cpr_use_drs_) { + OpmLog::warning("DRS_DISABLE", "Disabling DRS as matrix does not contain well contributions"); + } + parameters_.cpr_use_drs_ = false; if (parameters_.scale_linear_system_) { // also scale weights this->scaleEquationsAndVariables(weights_); @@ -288,23 +292,23 @@ protected: //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables //to be certain that correct matrix is used for preconditioning. Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel, - parallelInformation_ ); + parallelInformation_ ); assert( opA.comm() ); solve( opA, x, *rhs_, *(opA.comm()) ); } else { - typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false > Operator; - Operator opA(*matrix_, *matrix_, wellModel); + const WellModel& wellModel = simulator_.problem().wellModel(); + typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false > Operator; + Operator opA(*matrix_, *matrix_, wellModel); solve( opA, x, *rhs_ ); } - if (parameters_.scale_linear_system_) { + if (parameters_.scale_linear_system_) { scaleSolution(x); - } + } return converged_; - } @@ -404,13 +408,13 @@ protected: } - // 3x3 matrix block inversion was unstable at least 2.3 until and including - // 2.5.0. There may still be some issue with the 4x4 matrix block inversion - // we therefore still use the block inversion in OPM + // 3x3 matrix block inversion was unstable at least 2.3 until and including + // 2.5.0. There may still be some issue with the 4x4 matrix block inversion + // we therefore still use the block inversion in OPM typedef ParallelOverlappingILU0 >, - Vector, Vector> SeqPreconditioner; + Matrix::block_type::rows, + Matrix::block_type::cols> >, + Vector, Vector> SeqPreconditioner; template @@ -833,7 +837,7 @@ protected: std::vector>> overlapRowAndColumns_; FlowLinearSolverParameters parameters_; Vector weights_; - bool scale_variables_; + bool scale_variables_; }; // end ISTLSolver } // namespace Opm diff --git a/tests/test_blackoil_amg.cpp b/tests/test_blackoil_amg.cpp index d7ac54629..4fcbe7a9c 100644 --- a/tests/test_blackoil_amg.cpp +++ b/tests/test_blackoil_amg.cpp @@ -303,13 +303,8 @@ void runBlackoilAmgLaplace() smootherArgs.iterations = 1; Opm::CPRParameter param; - Vector weights(b.size()); - for (auto& elem : weights) { - elem = 1.0; - } - Opm::BlackoilAmg amg(param, - weights, + {}, fop, criterion, smootherArgs, comm);