From 40537f199914715ee5b03829c09ec213fca3c996 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 12 Mar 2019 13:55:11 +0100 Subject: [PATCH 01/17] Changes to make cpr work --- opm/autodiff/BlackoilAmg.hpp | 114 ++++--- opm/autodiff/CPRPreconditioner.hpp | 12 +- opm/autodiff/FlowLinearSolverParameters.hpp | 136 +++++--- opm/autodiff/ISTLSolverEbos.hpp | 331 +++++++++++++++++++- 4 files changed, 491 insertions(+), 102 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index 083dfbf94..be77d3dd3 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -21,6 +21,7 @@ #include #include +#include #include #include #include @@ -67,6 +68,11 @@ Dune::MatrixAdapter createOperator(const Dune::MatrixAdapter&, con return Dune::MatrixAdapter(matrix); } +template +std::unique_ptr< Dune::MatrixAdapter > createOperatorPtr(const Dune::MatrixAdapter&, const M& matrix, const T&) +{ + return std::make_unique< Dune::MatrixAdapter >(matrix); +} /** * \brief Creates an OverlappingSchwarzOperator as an operator. * @@ -89,30 +95,35 @@ Dune::OverlappingSchwarzOperator createOperator(const Dune::Overlapping //! \param comm The communication objecte describing the data distribution. //! \param pressureIndex The index of the pressure in the matrix block //! \retun A pair of the scaled matrix and the associated operator- -template +template std::tuple, Operator> -scaleMatrixQuasiImpes(const Operator& op, const Communication& comm, - std::size_t pressureIndex) +scaleMatrixDRS(const Operator& op, const Communication& comm, + std::size_t pressureIndex,const Vector& weights, const Opm::CPRParameter& param) { using Matrix = typename Operator::matrix_type; using Block = typename Matrix::block_type; + using BlockVector = typename Vector::block_type; std::unique_ptr matrix(new Matrix(op.getmat())); - - for ( auto& row : *matrix ) - { - for ( auto& block : row ) - { - for ( std::size_t i = 0; i < Block::rows; i++ ) - { - if ( i != pressureIndex ) - { - for(std::size_t j=0; j < Block::cols; j++) - { - block[pressureIndex][j] += block[i][j]; - } - } - } - } + if(param.cpr_use_drs_){ + const auto endi = matrix->end(); + for (auto i=matrix->begin(); i!=endi; ++i){ + const BlockVector& bw = weights[i.index()]; + const auto endj = (*i).end(); + 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][j] += block[i][j]; + } + } + block[pressureIndex] = bvec; + } + } + } } return std::make_tuple(std::move(matrix), createOperator(op, *matrix, comm)); } @@ -124,19 +135,20 @@ scaleMatrixQuasiImpes(const Operator& op, const Communication& comm, //! \param vector The vector to scale //! \param pressureIndex The index of the pressure in the matrix block template -void scaleVectorQuasiImpes(Vector& vector, std::size_t pressureIndex) +void scaleVectorDRS(Vector& vector, std::size_t pressureIndex, const Opm::CPRParameter& param, const Vector& weights) { using Block = typename Vector::block_type; - - for ( auto& block: vector) - { - for ( std::size_t i = 0; i < Block::dimension; i++ ) - { - if ( i != pressureIndex ) - { - block[pressureIndex] += block[i]; - } - } + 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] += block[i]; + } + block[pressureIndex] = val; + } } } @@ -400,10 +412,10 @@ private: } // Linear solver parameters const double tolerance = param_->cpr_solver_tol_; - const int maxit = param_->cpr_max_ell_iter_; + const int maxit = param_->cpr_max_iter_; const int verbosity = ( param_->cpr_solver_verbose_ && comm_.communicator().rank()==0 ) ? 1 : 0; - if ( param_->cpr_use_bicgstab_ ) + if ( param_->cpr_ell_solvetype_ == 0) { #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) Dune::BiCGSTABSolver solver(const_cast(op_), *sp, *prec, @@ -428,7 +440,7 @@ private: } #endif } - else + else if (param_->cpr_ell_solvetype_ == 1) { #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) Dune::CGSolver solver(const_cast(op_), *sp, *prec, @@ -453,6 +465,30 @@ private: } #endif } + else + { +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::LoopSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); + solver.apply(x,b,res); +#else + if ( !amg_ ) + { + Dune::LoopSolver solver(const_cast(op_), *sp, + reinterpret_cast(*prec), + tolerance, maxit, verbosity); + solver.apply(x,b,res); + } + else + { + Dune::LoopSolver solver(const_cast(op_), *sp, + reinterpret_cast(*prec), + tolerance, maxit, verbosity); + solver.apply(x,b,res); + } + +#endif + } #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) delete sp; @@ -873,7 +909,7 @@ private: * \brief An algebraic twolevel or multigrid approach for solving blackoil (supports CPR with and without AMG) * * This preconditioner first decouples the component used for coarsening using a simple scaling - * approach (e.g. Scheichl, Masson 2013,\see scaleMatrixQuasiImpes). Then it constructs the first + * approach (e.g. Scheichl, Masson 2013,\see scaleMatrixDRS). Then it constructs the first * coarse level system, either by simply extracting the coupling between the components at COMPONENT_INDEX * in the matrix blocks or by extracting them and applying aggregation to them directly. This coarse level * can be solved either by AMG or by ILU. The preconditioner is configured using CPRParameter. @@ -944,11 +980,13 @@ public: * \param comm The information about the parallelization. */ BlackoilAmg(const CPRParameter& param, + const typename TwoLevelMethod::FineDomainType& weights, const Operator& fineOperator, const Criterion& criterion, const SmootherArgs& smargs, const Communication& comm) : param_(param), - scaledMatrixOperator_(Detail::scaleMatrixQuasiImpes(fineOperator, comm, - COMPONENT_INDEX)), + weights_(weights), + scaledMatrixOperator_(Detail::scaleMatrixDRS(fineOperator, comm, + COMPONENT_INDEX, weights, param)), smoother_(Detail::constructSmoother(std::get<1>(scaledMatrixOperator_), smargs, comm)), levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_), @@ -973,16 +1011,18 @@ public: const typename TwoLevelMethod::FineRangeType& d) { auto scaledD = d; - Detail::scaleVectorQuasiImpes(scaledD, COMPONENT_INDEX); + Detail::scaleVectorDRS(scaledD, COMPONENT_INDEX, param_, weights_); twoLevelMethod_.apply(v, scaledD); } private: const CPRParameter& param_; + const typename TwoLevelMethod::FineDomainType& weights_; std::tuple, Operator> scaledMatrixOperator_; std::shared_ptr smoother_; LevelTransferPolicy levelTransferPolicy_; CoarseSolverPolicy coarseSolverPolicy_; TwoLevelMethod twoLevelMethod_; + //BlockVector weights_; }; namespace ISTLUtility diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index 3ddaf49ad..923044252 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -183,18 +183,22 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax, milu)); } -template < class C, class Op, class P, class S, std::size_t index > +template < class C, class Op, class P, class S, std::size_t index,class Vector> inline void createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::unique_ptr< BlackoilAmg >& amgPtr, - const CPRParameter& params) + const CPRParameter& params, + const Vector& weights) { using AMG = BlackoilAmg; + const int verbosity = ( params.cpr_solver_verbose_ && + comm.communicator().rank()==0 ) ? 1 : 0; + // TODO: revise choice of parameters int coarsenTarget=1200; using Criterion = C; Criterion criterion(15, coarsenTarget); - criterion.setDebugLevel( 0 ); // no debug information, 1 for printing hierarchy information + criterion.setDebugLevel( verbosity ); // no debug information, 1 for printing hierarchy information criterion.setDefaultValuesIsotropic(2); criterion.setNoPostSmoothSteps( 1 ); criterion.setNoPreSmoothSteps( 1 ); @@ -207,7 +211,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, smootherArgs.relaxationFactor = relax; setILUParameters(smootherArgs, params); - amgPtr.reset( new AMG( params, opA, criterion, smootherArgs, comm ) ); + amgPtr.reset( new AMG( params, weights, opA, criterion, smootherArgs, comm ) ); } template < class C, class Op, class P, class AMG > diff --git a/opm/autodiff/FlowLinearSolverParameters.hpp b/opm/autodiff/FlowLinearSolverParameters.hpp index 8f2e1fe98..0918eee15 100644 --- a/opm/autodiff/FlowLinearSolverParameters.hpp +++ b/opm/autodiff/FlowLinearSolverParameters.hpp @@ -59,6 +59,13 @@ NEW_PROP_TAG(UseAmg); NEW_PROP_TAG(UseCpr); NEW_PROP_TAG(LinearSolverBackend); NEW_PROP_TAG(PreconditionerAddWellContributions); +NEW_PROP_TAG(SystemStrategy); +NEW_PROP_TAG(ScaleLinearSystem); +NEW_PROP_TAG(CprSolverVerbose); +NEW_PROP_TAG(CprUseDrs); +NEW_PROP_TAG(CprMaxIter); +NEW_PROP_TAG(CprEllSolvetype); +NEW_PROP_TAG(CprReuseSetup); SET_SCALAR_PROP(FlowIstlSolverParams, LinearSolverReduction, 1e-2); SET_SCALAR_PROP(FlowIstlSolverParams, IluRelaxation, 0.9); @@ -76,6 +83,13 @@ SET_BOOL_PROP(FlowIstlSolverParams, UseAmg, false); SET_BOOL_PROP(FlowIstlSolverParams, UseCpr, false); SET_TYPE_PROP(FlowIstlSolverParams, LinearSolverBackend, Opm::ISTLSolverEbos); SET_BOOL_PROP(FlowIstlSolverParams, PreconditionerAddWellContributions, false); +SET_STRING_PROP(FlowIstlSolverParams, SystemStrategy, "original"); +SET_BOOL_PROP(FlowIstlSolverParams, ScaleLinearSystem, false); +SET_INT_PROP(FlowIstlSolverParams, CprSolverVerbose, 0); +SET_BOOL_PROP(FlowIstlSolverParams, CprUseDrs, false); +SET_INT_PROP(FlowIstlSolverParams, CprMaxIter, 20); +SET_INT_PROP(FlowIstlSolverParams, CprEllSolvetype, 0); +SET_INT_PROP(FlowIstlSolverParams, CprReuseSetup, 0); @@ -89,53 +103,60 @@ namespace Opm */ struct CPRParameter { - double cpr_relax_; - double cpr_solver_tol_; - int cpr_ilu_n_; + double cpr_relax_; + double cpr_solver_tol_; + int cpr_ilu_n_; MILU_VARIANT cpr_ilu_milu_; bool cpr_ilu_redblack_; bool cpr_ilu_reorder_sphere_; + bool cpr_use_drs_; int cpr_max_ell_iter_; + int cpr_max_iter_; + int cpr_ell_solvetype_; bool cpr_use_amg_; bool cpr_use_bicgstab_; - bool cpr_solver_verbose_; + int cpr_solver_verbose_; bool cpr_pressure_aggregation_; - + int cpr_reuse_setup_; CPRParameter() { reset(); } - CPRParameter( const ParameterGroup& param) - { - // reset values to default - reset(); + // CPRParameter( const ParameterGroup& param) + // { + // // reset values to default + // reset(); - cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_); - cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_); - cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_); - cpr_ilu_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_); - cpr_ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_); - cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_); - cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_); - cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_); - cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_); - cpr_pressure_aggregation_ = param.getDefault("cpr_pressure_aggregation", cpr_pressure_aggregation_); + // cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_); + // cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_); + // //cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_); + // //cpr_ilu_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_); + // //cpr_ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_); + // //cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_); + // //cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_); + // cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_); + // cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_); + // cpr_pressure_aggregation_ = param.getDefault("cpr_pressure_aggregation", cpr_pressure_aggregation_); - std::string milu("ILU"); - cpr_ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu)); - } + // //std::string milu("ILU"); + // //cpr_ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu)); + // } void reset() { - cpr_relax_ = 1.0; - cpr_solver_tol_ = 1e-2; + //cpr_relax_ = 1.0; + cpr_solver_tol_ = 1e-2; cpr_ilu_n_ = 0; cpr_ilu_milu_ = MILU_VARIANT::ILU; 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_use_amg_ = true; cpr_use_bicgstab_ = true; - cpr_solver_verbose_ = false; + cpr_solver_verbose_ = 0; cpr_pressure_aggregation_ = false; + cpr_reuse_setup_ = 0; } }; @@ -160,6 +181,8 @@ namespace Opm bool ignoreConvergenceFailure_; bool linear_solver_use_amg_; bool use_cpr_; + std::string system_strategy_; + bool scale_linear_system_; template void init() @@ -179,6 +202,14 @@ 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); + } template @@ -198,37 +229,44 @@ 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"); } FlowLinearSolverParameters() { reset(); } // read values from parameter class - FlowLinearSolverParameters( const ParameterGroup& param ) - : CPRParameter(param) - { - // set default parameters - reset(); + // FlowLinearSolverParameters( const ParameterGroup& param ) + // : CPRParameter(param) + // { + // // set default parameters + // reset(); - // read parameters (using previsouly set default values) - newton_use_gmres_ = param.getDefault("newton_use_gmres", newton_use_gmres_ ); - linear_solver_reduction_ = param.getDefault("linear_solver_reduction", linear_solver_reduction_ ); - linear_solver_maxiter_ = param.getDefault("linear_solver_maxiter", linear_solver_maxiter_); - linear_solver_restart_ = param.getDefault("linear_solver_restart", linear_solver_restart_); - linear_solver_verbosity_ = param.getDefault("linear_solver_verbosity", linear_solver_verbosity_); - require_full_sparsity_pattern_ = param.getDefault("require_full_sparsity_pattern", require_full_sparsity_pattern_); - ignoreConvergenceFailure_ = param.getDefault("linear_solver_ignoreconvergencefailure", ignoreConvergenceFailure_); - linear_solver_use_amg_ = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ ); - ilu_relaxation_ = param.getDefault("ilu_relaxation", ilu_relaxation_ ); - ilu_fillin_level_ = param.getDefault("ilu_fillin_level", ilu_fillin_level_ ); - ilu_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_); - ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_); - std::string milu("ILU"); - ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu)); + // // read parameters (using previsouly set default values) + // newton_use_gmres_ = param.getDefault("newton_use_gmres", newton_use_gmres_ ); + // linear_solver_reduction_ = param.getDefault("linear_solver_reduction", linear_solver_reduction_ ); + // linear_solver_maxiter_ = param.getDefault("linear_solver_maxiter", linear_solver_maxiter_); + // linear_solver_restart_ = param.getDefault("linear_solver_restart", linear_solver_restart_); + // linear_solver_verbosity_ = param.getDefault("linear_solver_verbosity", linear_solver_verbosity_); + // require_full_sparsity_pattern_ = param.getDefault("require_full_sparsity_pattern", require_full_sparsity_pattern_); + // ignoreConvergenceFailure_ = param.getDefault("linear_solver_ignoreconvergencefailure", ignoreConvergenceFailure_); + // linear_solver_use_amg_ = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ ); + // ilu_relaxation_ = param.getDefault("ilu_relaxation", ilu_relaxation_ ); + // ilu_fillin_level_ = param.getDefault("ilu_fillin_level", ilu_fillin_level_ ); + // ilu_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_); + // ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_); + // std::string milu("ILU"); + // ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu)); - // Check whether to use cpr approach - const std::string cprSolver = "cpr"; - use_cpr_ = ( param.getDefault("solver_approach", std::string()) == cprSolver ); - } + // // Check whether to use cpr approach + // const std::string cprSolver = "cpr"; + // use_cpr_ = ( param.getDefault("solver_approach", std::string()) == cprSolver ); + // } // set default values void reset() diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index e00fd6a16..73aa344fd 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -34,6 +34,7 @@ #include #include +//#include #include #include @@ -170,6 +171,7 @@ protected: template class ISTLSolverEbos { + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; @@ -178,10 +180,15 @@ protected: typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename SparseMatrixAdapter::IstlMatrix Matrix; + typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType; + typedef typename Vector::block_type BlockVector; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager; + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; enum { pressureIndex = Indices::pressureSwitchIdx }; static const int numEq = Indices::numEq; - - + public: typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType; @@ -207,12 +214,66 @@ protected: void eraseMatrix() { matrix_for_preconditioner_.reset(); } + void prepare(const SparseMatrixAdapter& M, Vector& b) { + matrix_.reset(new Matrix(M.istlMatrix())); + rhs_ = &b; + this->scaleSystem(); + } + void scaleSystem(){ + bool matrix_cont_added = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); + + if(matrix_cont_added){ + //Vector weights; + bool form_cpr = true; + if(parameters_.system_strategy_ == "quasiimpes"){ + weights_ = getQuasiImpesWeights(); + }else if(parameters_.system_strategy_ == "trueimpes"){ + weights_ = getStorageWeights(); + }else if(parameters_.system_strategy_ == "simple"){ + BlockVector bvec(1.0); + weights_ = getSimpleWeights(bvec); + }else if(parameters_.system_strategy_ == "original"){ + BlockVector bvec(0.0); + bvec[pressureIndex] = 1; + weights_ = getSimpleWeights(bvec); + }else{ + form_cpr = false; + } + // if(parameters_.linear_solver_verbosity_ > 1000) { + // std::ofstream filem("matrix_istl_pre.txt"); + // Dune::writeMatrixMarket(*matrix_, filem); + // std::ofstream fileb("rhs_istl_pre.txt"); + // Dune::writeMatrixMarket(*rhs_, fileb); + // std::ofstream filew("weights_istl.txt"); + // Dune::writeMatrixMarket(weights_, filew); + // } - void prepare(const SparseMatrixAdapter& M, const Vector& b) { - } + if(parameters_.scale_linear_system_){ + // also scale weights + this->scaleEquationsAndVariables(weights_); + } + if(form_cpr && not(parameters_.cpr_use_drs_)){ + scaleMatrixAndRhs(weights_); + } + + if(weights_.size() == 0){ + // if weights are not set cpr_use_drs_=false; + parameters_.cpr_use_drs_ = false; + } + + }else{ + if(parameters_.scale_linear_system_){ + // also scale weights + this->scaleEquationsAndVariables(weights_); + } + } + } + + + void setResidual(Vector& b) { - rhs_ = &b; + //rhs_ = &b; } void getResidual(Vector& b) const { @@ -220,7 +281,7 @@ protected: } void setMatrix(const SparseMatrixAdapter& M) { - matrix_ = &M.istlMatrix(); + //matrix_ = &M.istlMatrix(); } bool solve(Vector& x) { @@ -244,11 +305,50 @@ protected: 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_.linear_solver_verbosity_ > 5) && + // (iterations_ > parameters_.linear_solver_verbosity_)) { + // std::string dir = simulator_.problem().outputDir(); + // if (dir == ".") + // dir = ""; + // else if (!dir.empty() && dir.back() != '/') + // dir += "/"; + // namespace fs = boost::filesystem; + // fs::path output_dir(dir); + // fs::path subdir("reports"); + // output_dir = output_dir / subdir; + // if(!(fs::exists(output_dir))){ + // fs::create_directory(output_dir); + // } + // // Combine and return. + // std::ostringstream oss; + // oss << "prob_" << simulator_.episodeIndex() << "_"; + // oss << simulator_.time() << "_"; + // std::string output_file(oss.str()); + // fs::path full_path = output_dir / output_file; + // std::string prefix = full_path.string(); + // { + // std::string filename = prefix + "matrix_istl.txt"; + // std::ofstream filem(filename); + // Dune::writeMatrixMarket(*matrix_, filem); + // } + // { + // std::string filename = prefix + "rhs_istl.txt"; + // std::ofstream fileb(filename); + // Dune::writeMatrixMarket(*rhs_, fileb); + // } + // } + } + if(parameters_.scale_linear_system_){ + scaleSolution(x); + } + return converged_; @@ -411,7 +511,7 @@ protected: const MILU_VARIANT milu ) const { ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, - comm, amg, parameters_ ); + comm, amg, parameters_, weights_ ); } @@ -540,6 +640,7 @@ protected: protected: bool isParallel() const { + #if HAVE_MPI return parallelInformation_.type() == typeid(ParallelISTLInformation); #else @@ -574,17 +675,223 @@ protected: } } + // weights to make approxiate pressure equations + Vector getStorageWeights(){ + Vector weights(rhs_->size()); + BlockVector rhs(0.0); + rhs[pressureIndex] = 1.0; + int index = 0; + ElementContext elemCtx(simulator_); + const auto& vanguard = simulator_.vanguard(); + auto elemIt = vanguard.gridView().template begin(); + const auto& elemEndIt = vanguard.gridView().template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const Element& elem = *elemIt; + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + Dune::FieldVector storage; + unsigned threadId = ThreadManager::threadId(); + simulator_.model().localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,/*spaceIdx=*/0, /*timeIdx=*/0); + Scalar extrusionFactor = + elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor(); + Scalar scvVolume = + elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor; + Scalar storage_scale = scvVolume / elemCtx.simulator().timeStepSize(); + MatrixBlockType block; + int offset = 0; + double pressure_scale = 50e5; + for(int ii=0; ii< numEq; ++ii){ + for(int jj=0; jj< numEq; ++jj){ + //const auto& vec = storage[ii].derivative(jj); + block[ii][jj] = storage[ii].derivative(jj)/storage_scale; + if(jj==0){ + block[ii][jj] *=pressure_scale; + } + } + } + BlockVector bweights; + MatrixBlockType block_transpose = block.transpose(); + block_transpose.solve(bweights, rhs); + bweights /=1000; // given normal desnistyies this scales weights to about 1 + weights[index] = bweights; + ++index; + } + return weights; + } + + void scaleEquationsAndVariables(Vector& weights){ + // loop over primary variables + const auto& sol = simulator_.model().solution(0); + const auto endi = matrix_->end(); + int index = 0; + for (auto i=matrix_->begin(); i!=endi; ++i){ + const auto endj = (*i).end(); + BlockVector& brhs = (*rhs_)[i.index()]; + + for (auto j=(*i).begin(); j!=endj; ++j){ + MatrixBlockType& block = *j; + const auto& priVars = sol[i.index()]; + for ( std::size_t ii = 0; ii < block.rows; ii++ ){ + for(std::size_t jj=0; jj < block.cols; jj++){ + //double var_scale = getVarscale(jj, priVars.primaryVarsMeaning)) + double var_scale = simulator_.model().primaryVarWeight(i.index(),jj); + block[ii][jj] /=var_scale; + block[ii][jj] *= simulator_.model().eqWeight(i.index(), ii); + } + } + } + for(std::size_t ii=0; ii < brhs.size(); ii++){ + brhs[ii] *= simulator_.model().eqWeight(i.index(), ii); + } + if(weights_.size() == matrix_->N()){ + BlockVector& bw = weights[i.index()]; + for(std::size_t ii=0; ii < brhs.size(); ii++){ + bw[ii] /= simulator_.model().eqWeight(i.index(), ii); + } + double abs_max = + *std::max_element(bw.begin(), bw.end(), [](double a, double b){ return std::abs(a) < std::abs(b); } ); + bw /= abs_max; + } + + } + } + void scaleSolution(Vector& x){ + const auto& sol = simulator_.model().solution(0); + for(std::size_t i=0; i < x.size(); ++i){ + const auto& primVar = sol[i]; + auto& bx = x[i]; + for(std::size_t jj=0; jj < bx.size(); jj++){ + double var_scale = simulator_.model().primaryVarWeight(i,jj); + bx[jj] /= var_scale; + } + } + } + + Vector getQuasiImpesWeights(){ + Matrix& A = *matrix_; + Vector weights(rhs_->size()); + BlockVector rhs(0.0); + rhs[pressureIndex] = 1; + const auto endi = A.end(); + int index = 0; + for (auto i=A.begin(); i!=endi; ++i){ + const auto endj = (*i).end(); + MatrixBlockType diag_block(0.0); + for (auto j=(*i).begin(); j!=endj; ++j){ + if(i.index() == j.index()){ + diag_block = (*j); + break; + } + } + BlockVector bweights; + auto diag_block_transpose = diag_block.transpose(); + diag_block_transpose.solve(bweights, rhs); + double abs_max = + *std::max_element(bweights.begin(), bweights.end(), [](double a, double b){ return std::abs(a) < std::abs(b); } ); + bweights /= std::abs(abs_max); + weights[i.index()] = bweights; + } + return weights; + } + + Vector getSimpleWeights(const BlockVector& rhs){ + Vector weights(rhs_->size(),0); + for(auto& bw: weights){ + bw = rhs; + } + return weights; + } + + void scaleMatrixAndRhs(const Vector& weights){ + //static_assert(pressureIndex == 0, "Current support that pressure equation should be first"); + //using Matrix = typename Operator::matrix_type; + using Block = typename Matrix::block_type; + //Vector& rhs = *rhs_; + //Matrix& A = *matrix_; + //int index = 0; + //for ( auto& row : *matrix_ ){ + const auto endi = matrix_->end(); + for (auto i=matrix_->begin(); i!=endi; ++i){ + + //const auto& bweights = weights[row.index()]; + const BlockVector& bweights = weights[i.index()]; + BlockVector& brhs = (*rhs_)[i.index()]; + //++index; + //for ( auto& block : row ){ + const auto endj = (*i).end(); + for (auto j=(*i).begin(); j!=endj; ++j){ + // assume it is something on all rows + // the blew logic depend on pressureIndex=0 + Block& block = (*j); + for ( std::size_t ii = 0; ii < block.rows; ii++ ){ + if ( ii == 0 ){ + for(std::size_t jj=0; jj < block.cols; jj++){ + block[0][jj] *= bweights[ii];//*block[ii][jj]; + } + } else { + for(std::size_t jj=0; jj < block.cols; jj++){ + block[0][jj] += bweights[ii]*block[ii][jj]; + } + } + } + + } + for(std::size_t ii=0; ii < brhs.size(); ii++){ + if ( ii == 0 ){ + brhs[0] *= bweights[ii];//*brhs[ii]; + }else{ + brhs[0] += bweights[ii]*brhs[ii]; + } + } + } + } + + static void multBlocksInMatrix(Matrix& ebosJac,const MatrixBlockType& trans,bool left=true){ + const int n = ebosJac.N(); + //const int np = FluidSystem::numPhases; + for (int row_index = 0; row_index < n; ++row_index) { + auto& row = ebosJac[row_index]; + auto* dataptr = row.getptr(); + //auto* indexptr = row.getindexptr(); + for (int elem = 0; elem < row.N(); ++elem) { + auto& block = dataptr[elem]; + if(left){ + block = block.leftmultiply(trans); + }else{ + block = block.rightmultiply(trans); + } + } + } + } + + static void multBlocksVector(Vector& ebosResid_cp,const MatrixBlockType& leftTrans){ + for( auto& bvec: ebosResid_cp){ + auto bvec_new=bvec; + leftTrans.mv(bvec, bvec_new); + bvec=bvec_new; + } + } + static void scaleCPRSystem(Matrix& M_cp,Vector& b_cp,const MatrixBlockType& leftTrans){ + multBlocksInMatrix(M_cp, leftTrans, true); + multBlocksVector(b_cp, leftTrans); + } + + + const Simulator& simulator_; mutable int iterations_; mutable bool converged_; boost::any parallelInformation_; bool isIORank_; - const Matrix *matrix_; + + std::unique_ptr matrix_; Vector *rhs_; std::unique_ptr matrix_for_preconditioner_; std::vector>> overlapRowAndColumns_; FlowLinearSolverParameters parameters_; + Vector weights_; + bool scale_variables_; }; // end ISTLSolver } // namespace Opm From 3bcc80cc5b7c1d091ab7d15f1e85da6994658228 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 12 Mar 2019 15:13:38 +0100 Subject: [PATCH 02/17] changes to make test compile, however the tests tests matrixes functionality not used in flow --- tests/test_invert.cpp | 3 ++- tests/test_multmatrixtransposed.cpp | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/test_invert.cpp b/tests/test_invert.cpp index 76882bcf9..3ea71ba93 100644 --- a/tests/test_invert.cpp +++ b/tests/test_invert.cpp @@ -21,7 +21,8 @@ #define BOOST_TEST_MODULE InvertSpecializationTest #include -#include +#include + void checkIdentity(Dune::FieldMatrix M) { double diag = 0.0; diff --git a/tests/test_multmatrixtransposed.cpp b/tests/test_multmatrixtransposed.cpp index 19db6c646..0095a3c30 100644 --- a/tests/test_multmatrixtransposed.cpp +++ b/tests/test_multmatrixtransposed.cpp @@ -21,7 +21,7 @@ #define BOOST_TEST_MODULE MultMatrixTransposed #include -#include +#include using namespace Dune; using namespace Opm::Detail; From 13cd05fc1973843edce87eb1e4e27b8aa10fef5f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 14 Mar 2019 09:56:12 +0100 Subject: [PATCH 03/17] Add weights argument to test for BlackoilAmg. --- tests/test_blackoil_amg.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/tests/test_blackoil_amg.cpp b/tests/test_blackoil_amg.cpp index 259eb2f02..d7ac54629 100644 --- a/tests/test_blackoil_amg.cpp +++ b/tests/test_blackoil_amg.cpp @@ -303,7 +303,13 @@ 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); From 4eb44b4808546e059c763110e9620def96649733 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 14 Mar 2019 10:05:20 +0100 Subject: [PATCH 04/17] Remove unused code and fix indentation changes. --- opm/autodiff/FlowLinearSolverParameters.hpp | 56 ++------------------- 1 file changed, 4 insertions(+), 52 deletions(-) diff --git a/opm/autodiff/FlowLinearSolverParameters.hpp b/opm/autodiff/FlowLinearSolverParameters.hpp index 0918eee15..7ab940cf8 100644 --- a/opm/autodiff/FlowLinearSolverParameters.hpp +++ b/opm/autodiff/FlowLinearSolverParameters.hpp @@ -103,9 +103,9 @@ namespace Opm */ struct CPRParameter { - double cpr_relax_; - double cpr_solver_tol_; - int cpr_ilu_n_; + double cpr_relax_; + double cpr_solver_tol_; + int cpr_ilu_n_; MILU_VARIANT cpr_ilu_milu_; bool cpr_ilu_redblack_; bool cpr_ilu_reorder_sphere_; @@ -120,30 +120,9 @@ namespace Opm int cpr_reuse_setup_; CPRParameter() { reset(); } - // CPRParameter( const ParameterGroup& param) - // { - // // reset values to default - // reset(); - - // cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_); - // cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_); - // //cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_); - // //cpr_ilu_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_); - // //cpr_ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_); - // //cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_); - // //cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_); - // cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_); - // cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_); - // cpr_pressure_aggregation_ = param.getDefault("cpr_pressure_aggregation", cpr_pressure_aggregation_); - - // //std::string milu("ILU"); - // //cpr_ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu)); - // } - void reset() { - //cpr_relax_ = 1.0; - cpr_solver_tol_ = 1e-2; + cpr_solver_tol_ = 1e-2; cpr_ilu_n_ = 0; cpr_ilu_milu_ = MILU_VARIANT::ILU; cpr_ilu_redblack_ = false; @@ -240,33 +219,6 @@ namespace Opm } FlowLinearSolverParameters() { reset(); } - // read values from parameter class - // FlowLinearSolverParameters( const ParameterGroup& param ) - // : CPRParameter(param) - // { - // // set default parameters - // reset(); - - // // read parameters (using previsouly set default values) - // newton_use_gmres_ = param.getDefault("newton_use_gmres", newton_use_gmres_ ); - // linear_solver_reduction_ = param.getDefault("linear_solver_reduction", linear_solver_reduction_ ); - // linear_solver_maxiter_ = param.getDefault("linear_solver_maxiter", linear_solver_maxiter_); - // linear_solver_restart_ = param.getDefault("linear_solver_restart", linear_solver_restart_); - // linear_solver_verbosity_ = param.getDefault("linear_solver_verbosity", linear_solver_verbosity_); - // require_full_sparsity_pattern_ = param.getDefault("require_full_sparsity_pattern", require_full_sparsity_pattern_); - // ignoreConvergenceFailure_ = param.getDefault("linear_solver_ignoreconvergencefailure", ignoreConvergenceFailure_); - // linear_solver_use_amg_ = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ ); - // ilu_relaxation_ = param.getDefault("ilu_relaxation", ilu_relaxation_ ); - // ilu_fillin_level_ = param.getDefault("ilu_fillin_level", ilu_fillin_level_ ); - // ilu_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_); - // ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_); - // std::string milu("ILU"); - // ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu)); - - // // Check whether to use cpr approach - // const std::string cprSolver = "cpr"; - // use_cpr_ = ( param.getDefault("solver_approach", std::string()) == cprSolver ); - // } // set default values void reset() From 6ee5406a9fd8d51a9f2270b5de942d38d2b45ea7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 14 Mar 2019 10:37:46 +0100 Subject: [PATCH 05/17] Formatting and indentation fixes. --- opm/autodiff/BlackoilAmg.hpp | 70 +++++++++++++++--------------- opm/autodiff/CPRPreconditioner.hpp | 5 +-- 2 files changed, 38 insertions(+), 37 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index be77d3dd3..e643f9ad4 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -68,11 +68,18 @@ Dune::MatrixAdapter createOperator(const Dune::MatrixAdapter&, con return Dune::MatrixAdapter(matrix); } +/** + * \brief Creates a MatrixAdapter as an operator, storing it in a unique_ptr. + * + * The first argument is used to specify the return type using function overloading. + * \param matrix The matrix to wrap. + */ template std::unique_ptr< Dune::MatrixAdapter > createOperatorPtr(const Dune::MatrixAdapter&, const M& matrix, const T&) { return std::make_unique< Dune::MatrixAdapter >(matrix); } + /** * \brief Creates an OverlappingSchwarzOperator as an operator. * @@ -95,35 +102,32 @@ Dune::OverlappingSchwarzOperator createOperator(const Dune::Overlapping //! \param comm The communication objecte describing the data distribution. //! \param pressureIndex The index of the pressure in the matrix block //! \retun A pair of the scaled matrix and the associated operator- -template +template std::tuple, Operator> scaleMatrixDRS(const Operator& op, const Communication& comm, - std::size_t pressureIndex,const Vector& weights, const Opm::CPRParameter& param) + std::size_t pressureIndex, const Vector& weights, const Opm::CPRParameter& param) { using Matrix = typename Operator::matrix_type; using Block = typename Matrix::block_type; using BlockVector = typename Vector::block_type; std::unique_ptr matrix(new Matrix(op.getmat())); - if(param.cpr_use_drs_){ - const auto endi = matrix->end(); - for (auto i=matrix->begin(); i!=endi; ++i){ - const BlockVector& bw = weights[i.index()]; - const auto endj = (*i).end(); - 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][j] += block[i][j]; - } - } - block[pressureIndex] = bvec; - } - } - } + if (param.cpr_use_drs_) { + const auto endi = matrix->end(); + for (auto i = matrix->begin(); i != endi; ++i) { + const BlockVector& bw = weights[i.index()]; + const auto endj = (*i).end(); + 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; + } + } } return std::make_tuple(std::move(matrix), createOperator(op, *matrix, comm)); } @@ -138,17 +142,16 @@ template void scaleVectorDRS(Vector& vector, std::size_t pressureIndex, const Opm::CPRParameter& param, const Vector& weights) { 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] += block[i]; - } - block[pressureIndex] = val; - } + 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; + } } } @@ -1022,7 +1025,6 @@ private: LevelTransferPolicy levelTransferPolicy_; CoarseSolverPolicy coarseSolverPolicy_; TwoLevelMethod twoLevelMethod_; - //BlockVector weights_; }; namespace ISTLUtility diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index 923044252..81dbff274 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -183,7 +183,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax, milu)); } -template < class C, class Op, class P, class S, std::size_t index,class Vector> +template < class C, class Op, class P, class S, std::size_t index, class Vector> inline void createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::unique_ptr< BlackoilAmg >& amgPtr, @@ -191,8 +191,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, const Vector& weights) { using AMG = BlackoilAmg; - const int verbosity = ( params.cpr_solver_verbose_ && - comm.communicator().rank()==0 ) ? 1 : 0; + const int verbosity = ( params.cpr_solver_verbose_ && comm.communicator().rank() == 0 ) ? 1 : 0; // TODO: revise choice of parameters int coarsenTarget=1200; From 111feead1424c146d49ea0349622f319bb0547b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 14 Mar 2019 11:22:07 +0100 Subject: [PATCH 06/17] Formatting fixes and removing unused code. --- opm/autodiff/ISTLSolverEbos.hpp | 558 ++++++++++++++------------------ 1 file changed, 250 insertions(+), 308 deletions(-) diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index 73aa344fd..9c10b7f50 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -34,7 +34,6 @@ #include #include -//#include #include #include @@ -171,7 +170,7 @@ protected: template class ISTLSolverEbos { - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; @@ -180,15 +179,15 @@ protected: typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename SparseMatrixAdapter::IstlMatrix Matrix; - typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType; - typedef typename Vector::block_type BlockVector; - typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; - typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager; - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType; + typedef typename Vector::block_type BlockVector; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager; + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; enum { pressureIndex = Indices::pressureSwitchIdx }; static const int numEq = Indices::numEq; - + public: typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType; @@ -214,74 +213,63 @@ protected: void eraseMatrix() { matrix_for_preconditioner_.reset(); } - void prepare(const SparseMatrixAdapter& M, Vector& b) { - matrix_.reset(new Matrix(M.istlMatrix())); - rhs_ = &b; - this->scaleSystem(); - } - void scaleSystem(){ - bool matrix_cont_added = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); + + void prepare(const SparseMatrixAdapter& M, Vector& b) + { + matrix_.reset(new Matrix(M.istlMatrix())); + rhs_ = &b; + this->scaleSystem(); + } + + void scaleSystem() + { + const bool matrix_cont_added = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); - if(matrix_cont_added){ - //Vector weights; - bool form_cpr = true; - if(parameters_.system_strategy_ == "quasiimpes"){ - weights_ = getQuasiImpesWeights(); - }else if(parameters_.system_strategy_ == "trueimpes"){ - weights_ = getStorageWeights(); - }else if(parameters_.system_strategy_ == "simple"){ - BlockVector bvec(1.0); - weights_ = getSimpleWeights(bvec); - }else if(parameters_.system_strategy_ == "original"){ - BlockVector bvec(0.0); - bvec[pressureIndex] = 1; - weights_ = getSimpleWeights(bvec); - }else{ - form_cpr = false; - } - // if(parameters_.linear_solver_verbosity_ > 1000) { - // std::ofstream filem("matrix_istl_pre.txt"); - // Dune::writeMatrixMarket(*matrix_, filem); - // std::ofstream fileb("rhs_istl_pre.txt"); - // Dune::writeMatrixMarket(*rhs_, fileb); - // std::ofstream filew("weights_istl.txt"); - // Dune::writeMatrixMarket(weights_, filew); - // } + if (matrix_cont_added) { + bool form_cpr = true; + if (parameters_.system_strategy_ == "quasiimpes") { + weights_ = getQuasiImpesWeights(); + } else if (parameters_.system_strategy_ == "trueimpes") { + weights_ = getStorageWeights(); + } else if (parameters_.system_strategy_ == "simple") { + BlockVector bvec(1.0); + weights_ = getSimpleWeights(bvec); + } else if (parameters_.system_strategy_ == "original") { + BlockVector bvec(0.0); + bvec[pressureIndex] = 1; + weights_ = getSimpleWeights(bvec); + } else { + form_cpr = false; + } + if (parameters_.scale_linear_system_) { + // also scale weights + this->scaleEquationsAndVariables(weights_); + } + if (form_cpr && not(parameters_.cpr_use_drs_)) { + scaleMatrixAndRhs(weights_); + } + if (weights_.size() == 0) { + // if weights are not set cpr_use_drs_=false; + parameters_.cpr_use_drs_ = false; + } + } else { + if (parameters_.scale_linear_system_) { + // also scale weights + this->scaleEquationsAndVariables(weights_); + } + } + } - if(parameters_.scale_linear_system_){ - // also scale weights - this->scaleEquationsAndVariables(weights_); - } - if(form_cpr && not(parameters_.cpr_use_drs_)){ - scaleMatrixAndRhs(weights_); - } - - if(weights_.size() == 0){ - // if weights are not set cpr_use_drs_=false; - parameters_.cpr_use_drs_ = false; - } - - }else{ - if(parameters_.scale_linear_system_){ - // also scale weights - this->scaleEquationsAndVariables(weights_); - } - } - } - - - - - void setResidual(Vector& b) { - //rhs_ = &b; + void setResidual(Vector& /* b */) { + // rhs_ = &b; // Must be handled in prepare() instead. } void getResidual(Vector& b) const { b = *rhs_; } - void setMatrix(const SparseMatrixAdapter& M) { - //matrix_ = &M.istlMatrix(); + void setMatrix(const SparseMatrixAdapter& /* M */) { + // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead. } bool solve(Vector& x) { @@ -305,50 +293,15 @@ protected: solve( opA, x, *rhs_, *(opA.comm()) ); } else - - { - 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_.linear_solver_verbosity_ > 5) && - // (iterations_ > parameters_.linear_solver_verbosity_)) { - // std::string dir = simulator_.problem().outputDir(); - // if (dir == ".") - // dir = ""; - // else if (!dir.empty() && dir.back() != '/') - // dir += "/"; - // namespace fs = boost::filesystem; - // fs::path output_dir(dir); - // fs::path subdir("reports"); - // output_dir = output_dir / subdir; - // if(!(fs::exists(output_dir))){ - // fs::create_directory(output_dir); - // } - // // Combine and return. - // std::ostringstream oss; - // oss << "prob_" << simulator_.episodeIndex() << "_"; - // oss << simulator_.time() << "_"; - // std::string output_file(oss.str()); - // fs::path full_path = output_dir / output_file; - // std::string prefix = full_path.string(); - // { - // std::string filename = prefix + "matrix_istl.txt"; - // std::ofstream filem(filename); - // Dune::writeMatrixMarket(*matrix_, filem); - // } - // { - // std::string filename = prefix + "rhs_istl.txt"; - // std::ofstream fileb(filename); - // Dune::writeMatrixMarket(*rhs_, fileb); - // } - // } - } - if(parameters_.scale_linear_system_){ - scaleSolution(x); + } + + if (parameters_.scale_linear_system_) { + scaleSolution(x); } - return converged_; @@ -640,7 +593,6 @@ protected: protected: bool isParallel() const { - #if HAVE_MPI return parallelInformation_.type() == typeid(ParallelISTLInformation); #else @@ -675,208 +627,198 @@ protected: } } - // weights to make approxiate pressure equations - Vector getStorageWeights(){ - Vector weights(rhs_->size()); - BlockVector rhs(0.0); - rhs[pressureIndex] = 1.0; - int index = 0; - ElementContext elemCtx(simulator_); - const auto& vanguard = simulator_.vanguard(); - auto elemIt = vanguard.gridView().template begin(); - const auto& elemEndIt = vanguard.gridView().template end(); - for (; elemIt != elemEndIt; ++elemIt) { - const Element& elem = *elemIt; - elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - Dune::FieldVector storage; - unsigned threadId = ThreadManager::threadId(); - simulator_.model().localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,/*spaceIdx=*/0, /*timeIdx=*/0); - Scalar extrusionFactor = - elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor(); - Scalar scvVolume = - elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor; - Scalar storage_scale = scvVolume / elemCtx.simulator().timeStepSize(); - MatrixBlockType block; - int offset = 0; - double pressure_scale = 50e5; - for(int ii=0; ii< numEq; ++ii){ - for(int jj=0; jj< numEq; ++jj){ - //const auto& vec = storage[ii].derivative(jj); - block[ii][jj] = storage[ii].derivative(jj)/storage_scale; - if(jj==0){ - block[ii][jj] *=pressure_scale; - } - } - } - BlockVector bweights; - MatrixBlockType block_transpose = block.transpose(); - block_transpose.solve(bweights, rhs); - bweights /=1000; // given normal desnistyies this scales weights to about 1 - weights[index] = bweights; - ++index; - } - return weights; - } + // Weights to make approximate pressure equations. + Vector getStorageWeights() const + { + Vector weights(rhs_->size()); + BlockVector rhs(0.0); + rhs[pressureIndex] = 1.0; + int index = 0; + ElementContext elemCtx(simulator_); + const auto& vanguard = simulator_.vanguard(); + auto elemIt = vanguard.gridView().template begin(); + const auto& elemEndIt = vanguard.gridView().template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const Element& elem = *elemIt; + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + Dune::FieldVector storage; + unsigned threadId = ThreadManager::threadId(); + simulator_.model().localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,/*spaceIdx=*/0, /*timeIdx=*/0); + Scalar extrusionFactor = elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor(); + Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor; + Scalar storage_scale = scvVolume / elemCtx.simulator().timeStepSize(); + MatrixBlockType block; + int offset = 0; + double pressure_scale = 50e5; + for (int ii = 0; ii < numEq; ++ii) { + for (int jj = 0; jj < numEq; ++jj) { + block[ii][jj] = storage[ii].derivative(jj)/storage_scale; + if (jj == 0) { + block[ii][jj] *= pressure_scale; + } + } + } + BlockVector bweights; + MatrixBlockType block_transpose = block.transpose(); + block_transpose.solve(bweights, rhs); + bweights /= 1000.0; // given normal densities this scales weights to about 1. + weights[index] = bweights; + ++index; + } + return weights; + } - void scaleEquationsAndVariables(Vector& weights){ - // loop over primary variables - const auto& sol = simulator_.model().solution(0); - const auto endi = matrix_->end(); - int index = 0; - for (auto i=matrix_->begin(); i!=endi; ++i){ - const auto endj = (*i).end(); - BlockVector& brhs = (*rhs_)[i.index()]; - - for (auto j=(*i).begin(); j!=endj; ++j){ - MatrixBlockType& block = *j; - const auto& priVars = sol[i.index()]; - for ( std::size_t ii = 0; ii < block.rows; ii++ ){ - for(std::size_t jj=0; jj < block.cols; jj++){ - //double var_scale = getVarscale(jj, priVars.primaryVarsMeaning)) - double var_scale = simulator_.model().primaryVarWeight(i.index(),jj); - block[ii][jj] /=var_scale; - block[ii][jj] *= simulator_.model().eqWeight(i.index(), ii); - } - } - } - for(std::size_t ii=0; ii < brhs.size(); ii++){ - brhs[ii] *= simulator_.model().eqWeight(i.index(), ii); - } - if(weights_.size() == matrix_->N()){ - BlockVector& bw = weights[i.index()]; - for(std::size_t ii=0; ii < brhs.size(); ii++){ - bw[ii] /= simulator_.model().eqWeight(i.index(), ii); - } - double abs_max = - *std::max_element(bw.begin(), bw.end(), [](double a, double b){ return std::abs(a) < std::abs(b); } ); - bw /= abs_max; - } - - } - } - void scaleSolution(Vector& x){ - const auto& sol = simulator_.model().solution(0); - for(std::size_t i=0; i < x.size(); ++i){ - const auto& primVar = sol[i]; - auto& bx = x[i]; - for(std::size_t jj=0; jj < bx.size(); jj++){ - double var_scale = simulator_.model().primaryVarWeight(i,jj); - bx[jj] /= var_scale; - } - } - } - - Vector getQuasiImpesWeights(){ - Matrix& A = *matrix_; - Vector weights(rhs_->size()); - BlockVector rhs(0.0); - rhs[pressureIndex] = 1; - const auto endi = A.end(); - int index = 0; - for (auto i=A.begin(); i!=endi; ++i){ - const auto endj = (*i).end(); - MatrixBlockType diag_block(0.0); - for (auto j=(*i).begin(); j!=endj; ++j){ - if(i.index() == j.index()){ - diag_block = (*j); - break; - } - } - BlockVector bweights; - auto diag_block_transpose = diag_block.transpose(); - diag_block_transpose.solve(bweights, rhs); - double abs_max = - *std::max_element(bweights.begin(), bweights.end(), [](double a, double b){ return std::abs(a) < std::abs(b); } ); - bweights /= std::abs(abs_max); - weights[i.index()] = bweights; - } - return weights; - } - - Vector getSimpleWeights(const BlockVector& rhs){ - Vector weights(rhs_->size(),0); - for(auto& bw: weights){ - bw = rhs; - } - return weights; - } + void scaleEquationsAndVariables(Vector& weights) + { + // loop over primary variables + const auto& sol = simulator_.model().solution(0); + const auto endi = matrix_->end(); + int index = 0; + for (auto i = matrix_->begin(); i != endi; ++i) { + const auto endj = (*i).end(); + BlockVector& brhs = (*rhs_)[i.index()]; + for (auto j = (*i).begin(); j != endj; ++j) { + MatrixBlockType& block = *j; + const auto& priVars = sol[i.index()]; + for (std::size_t ii = 0; ii < block.rows; ii++ ) { + for (std::size_t jj = 0; jj < block.cols; jj++) { + double var_scale = simulator_.model().primaryVarWeight(i.index(),jj); + block[ii][jj] /= var_scale; + block[ii][jj] *= simulator_.model().eqWeight(i.index(), ii); + } + } + } + for (std::size_t ii = 0; ii < brhs.size(); ii++) { + brhs[ii] *= simulator_.model().eqWeight(i.index(), ii); + } + if (weights_.size() == matrix_->N()) { + BlockVector& bw = weights[i.index()]; + for (std::size_t ii = 0; ii < brhs.size(); ii++) { + bw[ii] /= simulator_.model().eqWeight(i.index(), ii); + } + double abs_max = + *std::max_element(bw.begin(), bw.end(), [](double a, double b){ return std::abs(a) < std::abs(b); } ); + bw /= abs_max; + } + } + } - void scaleMatrixAndRhs(const Vector& weights){ - //static_assert(pressureIndex == 0, "Current support that pressure equation should be first"); - //using Matrix = typename Operator::matrix_type; - using Block = typename Matrix::block_type; - //Vector& rhs = *rhs_; - //Matrix& A = *matrix_; - //int index = 0; - //for ( auto& row : *matrix_ ){ - const auto endi = matrix_->end(); - for (auto i=matrix_->begin(); i!=endi; ++i){ - - //const auto& bweights = weights[row.index()]; - const BlockVector& bweights = weights[i.index()]; - BlockVector& brhs = (*rhs_)[i.index()]; - //++index; - //for ( auto& block : row ){ - const auto endj = (*i).end(); - for (auto j=(*i).begin(); j!=endj; ++j){ - // assume it is something on all rows - // the blew logic depend on pressureIndex=0 - Block& block = (*j); - for ( std::size_t ii = 0; ii < block.rows; ii++ ){ - if ( ii == 0 ){ - for(std::size_t jj=0; jj < block.cols; jj++){ - block[0][jj] *= bweights[ii];//*block[ii][jj]; - } - } else { - for(std::size_t jj=0; jj < block.cols; jj++){ - block[0][jj] += bweights[ii]*block[ii][jj]; - } - } - } - - } - for(std::size_t ii=0; ii < brhs.size(); ii++){ - if ( ii == 0 ){ - brhs[0] *= bweights[ii];//*brhs[ii]; - }else{ - brhs[0] += bweights[ii]*brhs[ii]; - } - } - } - } - - static void multBlocksInMatrix(Matrix& ebosJac,const MatrixBlockType& trans,bool left=true){ - const int n = ebosJac.N(); - //const int np = FluidSystem::numPhases; - for (int row_index = 0; row_index < n; ++row_index) { - auto& row = ebosJac[row_index]; - auto* dataptr = row.getptr(); - //auto* indexptr = row.getindexptr(); - for (int elem = 0; elem < row.N(); ++elem) { - auto& block = dataptr[elem]; - if(left){ - block = block.leftmultiply(trans); - }else{ - block = block.rightmultiply(trans); - } - } - } - } - - static void multBlocksVector(Vector& ebosResid_cp,const MatrixBlockType& leftTrans){ - for( auto& bvec: ebosResid_cp){ - auto bvec_new=bvec; - leftTrans.mv(bvec, bvec_new); - bvec=bvec_new; - } - } - static void scaleCPRSystem(Matrix& M_cp,Vector& b_cp,const MatrixBlockType& leftTrans){ - multBlocksInMatrix(M_cp, leftTrans, true); - multBlocksVector(b_cp, leftTrans); - } + void scaleSolution(Vector& x) + { + const auto& sol = simulator_.model().solution(0); + for (std::size_t i = 0; i < x.size(); ++i) { + const auto& primVar = sol[i]; + auto& bx = x[i]; + for (std::size_t jj = 0; jj < bx.size(); jj++) { + double var_scale = simulator_.model().primaryVarWeight(i,jj); + bx[jj] /= var_scale; + } + } + } + Vector getQuasiImpesWeights() + { + Matrix& A = *matrix_; + Vector weights(rhs_->size()); + BlockVector rhs(0.0); + rhs[pressureIndex] = 1; + const auto endi = A.end(); + int index = 0; + for (auto i = A.begin(); i!=endi; ++i) { + const auto endj = (*i).end(); + MatrixBlockType diag_block(0.0); + for (auto j=(*i).begin(); j!=endj; ++j) { + if (i.index() == j.index()) { + diag_block = (*j); + break; + } + } + BlockVector bweights; + auto diag_block_transpose = diag_block.transpose(); + diag_block_transpose.solve(bweights, rhs); + double abs_max = + *std::max_element(bweights.begin(), bweights.end(), [](double a, double b){ return std::abs(a) < std::abs(b); } ); + bweights /= std::abs(abs_max); + weights[i.index()] = bweights; + } + return weights; + } + Vector getSimpleWeights(const BlockVector& rhs) + { + Vector weights(rhs_->size(), 0); + for (auto& bw : weights) { + bw = rhs; + } + return weights; + } + + void scaleMatrixAndRhs(const Vector& weights) + { + using Block = typename Matrix::block_type; + const auto endi = matrix_->end(); + for (auto i = matrix_->begin(); i !=endi; ++i) { + const BlockVector& bweights = weights[i.index()]; + BlockVector& brhs = (*rhs_)[i.index()]; + const auto endj = (*i).end(); + for (auto j = (*i).begin(); j != endj; ++j) { + // assume it is something on all rows + // the blew logic depend on pressureIndex=0 + Block& block = (*j); + for ( std::size_t ii = 0; ii < block.rows; ii++ ) { + if ( ii == 0 ) { + for (std::size_t jj = 0; jj < block.cols; jj++) { + block[0][jj] *= bweights[ii]; + } + } else { + for (std::size_t jj = 0; jj < block.cols; jj++) { + block[0][jj] += bweights[ii]*block[ii][jj]; + } + } + } + } + for (std::size_t ii = 0; ii < brhs.size(); ii++) { + if ( ii == 0 ){ + brhs[0] *= bweights[ii]; + } else { + brhs[0] += bweights[ii]*brhs[ii]; + } + } + } + } + + static void multBlocksInMatrix(Matrix& ebosJac, const MatrixBlockType& trans, const bool left = true) + { + const int n = ebosJac.N(); + for (int row_index = 0; row_index < n; ++row_index) { + auto& row = ebosJac[row_index]; + auto* dataptr = row.getptr(); + for (int elem = 0; elem < row.N(); ++elem) { + auto& block = dataptr[elem]; + if (left) { + block = block.leftmultiply(trans); + } else { + block = block.rightmultiply(trans); + } + } + } + } + + static void multBlocksVector(Vector& ebosResid_cp, const MatrixBlockType& leftTrans) + { + for (auto& bvec : ebosResid_cp) { + auto bvec_new = bvec; + leftTrans.mv(bvec, bvec_new); + bvec = bvec_new; + } + } + + static void scaleCPRSystem(Matrix& M_cp, Vector& b_cp, const MatrixBlockType& leftTrans) + { + multBlocksInMatrix(M_cp, leftTrans, true); + multBlocksVector(b_cp, leftTrans); + } const Simulator& simulator_; mutable int iterations_; From e05e5e4f0e1337e8e2c25da15940f72ca712dc08 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 14 Mar 2019 09:04:51 +0100 Subject: [PATCH 07/17] Fix indenting --- opm/autodiff/BlackoilAmg.hpp | 26 +++++++------- opm/autodiff/CPRPreconditioner.hpp | 2 +- opm/autodiff/FlowLinearSolverParameters.hpp | 38 ++++++++++----------- opm/autodiff/ISTLSolverEbos.hpp | 26 +++++++------- 4 files changed, 44 insertions(+), 48 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index e643f9ad4..c5fda483c 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. * @@ -468,8 +467,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 +488,9 @@ private: tolerance, maxit, verbosity); solver.apply(x,b,res); } - -#endif - } + +#endif + } #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) delete sp; @@ -983,20 +982,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..2802f047e 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -288,23 +288,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 +404,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 +833,7 @@ protected: std::vector>> overlapRowAndColumns_; FlowLinearSolverParameters parameters_; Vector weights_; - bool scale_variables_; + bool scale_variables_; }; // end ISTLSolver } // namespace Opm From 41ff0da2b4ddb98d631bc4e04ece352bf19c0f86 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 14 Mar 2019 09:32:32 +0100 Subject: [PATCH 08/17] fixed: drs without matrix-add-well-contributions --- opm/autodiff/ISTLSolverEbos.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index 2802f047e..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_); From e12b959defcb4cacfdea44588b4a95cae98d05fa Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 14 Mar 2019 09:42:34 +0100 Subject: [PATCH 09/17] fix test_blackoil_amg --- tests/test_blackoil_amg.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) 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); From f98be9d43bb0d68c95430afe973b716bbcefa413 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 14 Mar 2019 11:25:45 +0100 Subject: [PATCH 10/17] simplify: use matrix-vector product in DRS for matrix --- opm/autodiff/BlackoilAmg.hpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index c5fda483c..9ce6de845 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -118,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); } } } From e1db2d2344eed0c6b88bda168479150d14c2d0e9 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 14 Mar 2019 11:27:21 +0100 Subject: [PATCH 11/17] simplify: use dotproduct in DRS for vector --- opm/autodiff/BlackoilAmg.hpp | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index 9ce6de845..ee18a087c 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -139,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); } } } From bedd88bc771988af6e51a659a18e2256b4b5a817 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 14 Mar 2019 15:21:34 +0100 Subject: [PATCH 12/17] remove unused variables --- opm/autodiff/BlackoilAmg.hpp | 1 - opm/autodiff/ISTLSolverEbos.hpp | 7 ------- 2 files changed, 8 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index ee18a087c..0e9288060 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -116,7 +116,6 @@ scaleMatrixDRS(const Operator& op, const Communication& comm, const BlockVector& bw = weights[i.index()]; const auto endj = (*i).end(); for (auto j = (*i).begin(); j != endj; ++j) { - BlockVector bvec(0.0); Block& block = *j; BlockVector& bvec = block[pressureIndex]; // should introduce limits which also change the weights diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index d59a02b88..0fc1c9015 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -653,7 +653,6 @@ protected: Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor; Scalar storage_scale = scvVolume / elemCtx.simulator().timeStepSize(); MatrixBlockType block; - int offset = 0; double pressure_scale = 50e5; for (int ii = 0; ii < numEq; ++ii) { for (int jj = 0; jj < numEq; ++jj) { @@ -676,15 +675,12 @@ protected: void scaleEquationsAndVariables(Vector& weights) { // loop over primary variables - const auto& sol = simulator_.model().solution(0); const auto endi = matrix_->end(); - int index = 0; for (auto i = matrix_->begin(); i != endi; ++i) { const auto endj = (*i).end(); BlockVector& brhs = (*rhs_)[i.index()]; for (auto j = (*i).begin(); j != endj; ++j) { MatrixBlockType& block = *j; - const auto& priVars = sol[i.index()]; for (std::size_t ii = 0; ii < block.rows; ii++ ) { for (std::size_t jj = 0; jj < block.cols; jj++) { double var_scale = simulator_.model().primaryVarWeight(i.index(),jj); @@ -710,9 +706,7 @@ protected: void scaleSolution(Vector& x) { - const auto& sol = simulator_.model().solution(0); for (std::size_t i = 0; i < x.size(); ++i) { - const auto& primVar = sol[i]; auto& bx = x[i]; for (std::size_t jj = 0; jj < bx.size(); jj++) { double var_scale = simulator_.model().primaryVarWeight(i,jj); @@ -728,7 +722,6 @@ protected: BlockVector rhs(0.0); rhs[pressureIndex] = 1; const auto endi = A.end(); - int index = 0; for (auto i = A.begin(); i!=endi; ++i) { const auto endj = (*i).end(); MatrixBlockType diag_block(0.0); From 31f4c0193883fc6bc282bfadce11ae6cf2022e2f Mon Sep 17 00:00:00 2001 From: hnil Date: Fri, 15 Mar 2019 11:27:38 +0100 Subject: [PATCH 13/17] Changes to make CPR able to use general indexing. --- opm/autodiff/BlackoilAmg.hpp | 67 ++++++++++++++++++++---------- opm/autodiff/CPRPreconditioner.hpp | 18 +++++--- opm/autodiff/ISTLSolverEbos.hpp | 48 +++++++++------------ tests/test_blackoil_amg.cpp | 2 +- 4 files changed, 78 insertions(+), 57 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index 0e9288060..88b3453a8 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -42,6 +42,27 @@ class UnSymmetricCriterion; } } +namespace Opm +{ + namespace Amg + { + template + class Element + { + public: + enum { /* @brief We preserve the sign.*/ + is_sign_preserving = true + }; + + template + typename M::field_type operator()(const M& m) const + { + return m[Row][Column]; + } + }; + } // namespace Amg +} // namespace Opm + namespace Dune { @@ -99,12 +120,12 @@ Dune::OverlappingSchwarzOperator createOperator(const Dune::Overlapping //! Sedimentary Basin Simulations, 2003. //! \param op The operator that stems from the discretization. //! \param comm The communication objecte describing the data distribution. -//! \param pressureIndex The index of the pressure in the matrix block +//! \param pressureEqnIndex The index of the pressure in the matrix block //! \retun A pair of the scaled matrix and the associated operator- template std::tuple, Operator> scaleMatrixDRS(const Operator& op, const Communication& comm, - std::size_t pressureIndex, const Vector& weights, const Opm::CPRParameter& param) + std::size_t pressureEqnIndex, const Vector& weights, const Opm::CPRParameter& param) { using Matrix = typename Operator::matrix_type; using Block = typename Matrix::block_type; @@ -117,7 +138,7 @@ scaleMatrixDRS(const Operator& op, const Communication& comm, const auto endj = (*i).end(); for (auto j = (*i).begin(); j != endj; ++j) { Block& block = *j; - BlockVector& bvec = block[pressureIndex]; + BlockVector& bvec = block[pressureEqnIndex]; // should introduce limits which also change the weights block.mtv(bw, bvec); } @@ -131,16 +152,16 @@ scaleMatrixDRS(const Operator& op, const Communication& comm, //! See section 3.2.3 of Scheichl, Masson: Decoupling and Block Preconditioning for //! Sedimentary Basin Simulations, 2003. //! \param vector The vector to scale -//! \param pressureIndex The index of the pressure in the matrix block +//! \param pressureEqnIndex The index of the pressure in the matrix block template -void scaleVectorDRS(Vector& vector, std::size_t pressureIndex, const Opm::CPRParameter& param, const Vector& weights) +void scaleVectorDRS(Vector& vector, std::size_t pressureEqnIndex, const Opm::CPRParameter& param, const Vector& weights) { using Block = typename Vector::block_type; if (param.cpr_use_drs_) { for (std::size_t j = 0; j < vector.size(); ++j) { Block& block = vector[j]; const Block& bw = weights[j]; - block[pressureIndex] = bw.dot(block); + block[pressureEqnIndex] = bw.dot(block); } } } @@ -289,23 +310,23 @@ struct ScalarType::value>, Dune::Amg::FirstDiagonal> >; }; -template +template struct OneComponentCriterionType {}; -template -struct OneComponentCriterionType,N> >,COMPONENT_INDEX> +template +struct OneComponentCriterionType,N> >, COMPONENT_INDEX, VARIABLE_INDEX> { - using value = Dune::Amg::CoarsenCriterion, Dune::Amg::Diagonal > >; + using value = Dune::Amg::CoarsenCriterion, Opm::Amg::Element > >; }; -template -struct OneComponentCriterionType,N> >,COMPONENT_INDEX> +template +struct OneComponentCriterionType,N> >, COMPONENT_INDEX, VARIABLE_INDEX> { - using value = Dune::Amg::CoarsenCriterion, Dune::Amg::Diagonal > >; + using value = Dune::Amg::CoarsenCriterion, Opm::Amg::Element > >; }; -template +template class OneComponentAggregationLevelTransferPolicy; @@ -674,7 +695,7 @@ void buildCoarseSparseMatrix(M& coarseMatrix, G& fineGraph, const V& visitedMap, * @tparam Criterion The criterion that describes the aggregation procedure. * @tparam Communication The class that describes the communication pattern. */ -template +template class OneComponentAggregationLevelTransferPolicy : public Dune::Amg::LevelTransferPolicy::value> { @@ -785,7 +806,7 @@ public: for ( auto col = row.begin(), cend = row.end(); col != cend; ++col, ++coarseCol ) { assert( col.index() == coarseCol.index() ); - *coarseCol = (*col)[COMPONENT_INDEX][COMPONENT_INDEX]; + *coarseCol = (*col)[COMPONENT_INDEX][VARIABLE_INDEX]; } ++coarseRow; } @@ -815,7 +836,7 @@ public: const auto& j = (*aggregatesMap_)[entry.index()]; if ( j != AggregatesMap::ISOLATED ) { - (*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][COMPONENT_INDEX]; + (*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][VARIABLE_INDEX]; } } } @@ -911,9 +932,10 @@ private: * \tparam C The type of coarsening criterion to use. * \tparam P The type of the class describing the parallelization. * \tparam COMPONENT_INDEX The index of the component to use for coarsening (usually the pressure). + * \tparam VARIABLE_INDEX The index of the variable to use for coarsening (usually the pressure). */ template + typename P, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX> class BlackoilAmg : public Dune::Preconditioner { @@ -934,13 +956,14 @@ protected: using CoarseOperator = typename Detail::ScalarType::value; using CoarseSmoother = typename Detail::ScalarType::value; using FineCriterion = - typename Detail::OneComponentCriterionType::value; + typename Detail::OneComponentCriterionType::value; using CoarseCriterion = typename Detail::ScalarType::value; using LevelTransferPolicy = OneComponentAggregationLevelTransferPolicy; + COMPONENT_INDEX, + VARIABLE_INDEX>; using CoarseSolverPolicy = Detail::OneStepAMGCoarseSolverPolicy +template struct BlackoilAmgSelector { using Criterion = C; @@ -1036,7 +1059,7 @@ struct BlackoilAmgSelector using ParallelInformation = typename Selector::ParallelInformation; using Operator = typename Selector::Operator; using Smoother = typename Selector::EllipticPreconditioner; - using AMG = BlackoilAmg; + using AMG = BlackoilAmg; }; } // end namespace ISTLUtility } // end namespace Opm diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index 1e46558f0..5ea90e94a 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -55,9 +55,15 @@ namespace Opm { template + typename P, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX> class BlackoilAmg; +namespace Amg +{ + template + class Element; +} + namespace ISTLUtility { @@ -183,14 +189,14 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax, milu)); } -template < class C, class Op, class P, class S, std::size_t index, class Vector> +template < class C, class Op, class P, class S, std::size_t PressureEqnIndex, std::size_t PressureVarIndex, class Vector> inline void createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, - std::unique_ptr< BlackoilAmg >& amgPtr, + std::unique_ptr< BlackoilAmg >& amgPtr, const CPRParameter& params, const Vector& weights) { - using AMG = BlackoilAmg; + using AMG = BlackoilAmg; const int verbosity = ( params.cpr_solver_verbose_ && comm.communicator().rank() == 0 ) ? 1 : 0; // TODO: revise choice of parameters @@ -242,7 +248,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const MILU_VARIANT m /// \param relax The relaxation parameter for ILU0. /// \param comm The object describing the parallelization information and communication. // \param amgPtr The unique_ptr to be filled (return) -template < int pressureIndex=0, class Op, class P, class AMG > +template < int PressureEqnIndex, int PressureVarIndex, class Op, class P, class AMG > inline void createAMGPreconditionerPointer( Op& opA, const double relax, const MILU_VARIANT milu, const P& comm, std::unique_ptr< AMG >& amgPtr ) { @@ -250,7 +256,7 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const MILU_VARIANT typedef typename Op::matrix_type M; // The coupling metric used in the AMG - typedef Dune::Amg::Diagonal CouplingMetric; + typedef Opm::Amg::Element CouplingMetric; // The coupling criterion used in the AMG typedef Dune::Amg::SymmetricCriterion CritBase; diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index 0fc1c9015..3440956d3 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -1,5 +1,6 @@ /* Copyright 2016 IRIS AS + Copyright 2019 Equinor ASA This file is part of the Open Porous Media project (OPM). @@ -31,6 +32,7 @@ #include #include #include +#include #include #include @@ -162,11 +164,6 @@ protected: /// solving the reduced system (after eliminating well variables) /// as a block-structured matrix (one block for all cell variables) for a fixed /// number of cell variables np . - /// \tparam MatrixBlockType The type of the matrix block used. - /// \tparam VectorBlockType The type of the vector block used. - /// \tparam pressureIndex The index of the pressure component in the vector - /// vector block. It is used to guide the AMG coarsening. - /// Default is zero. template class ISTLSolverEbos { @@ -185,7 +182,9 @@ protected: typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager; typedef typename GridView::template Codim<0>::Entity Element; typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; - enum { pressureIndex = Indices::pressureSwitchIdx }; + // Due to miscibility oil <-> gas the water eqn is the one we can replace with a pressure equation. + enum { pressureEqnIndex = BlackOilDefaultIndexTraits::waterCompIdx }; + enum { pressureVarIndex = Indices::pressureSwitchIdx }; static const int numEq = Indices::numEq; public: @@ -236,7 +235,7 @@ protected: weights_ = getSimpleWeights(bvec); } else if (parameters_.system_strategy_ == "original") { BlockVector bvec(0.0); - bvec[pressureIndex] = 1; + bvec[pressureEqnIndex] = 1; weights_ = getSimpleWeights(bvec); } else { form_cpr = false; @@ -370,11 +369,11 @@ protected: if ( parameters_.use_cpr_ ) { using Matrix = typename MatrixOperator::matrix_type; - using CouplingMetric = Dune::Amg::Diagonal; + using CouplingMetric = Opm::Amg::Element; using CritBase = Dune::Amg::SymmetricCriterion; using Criterion = Dune::Amg::CoarsenCriterion; using AMG = typename ISTLUtility - ::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureIndex >::AMG; + ::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureEqnIndex, pressureVarIndex >::AMG; std::unique_ptr< AMG > amg; // Construct preconditioner. @@ -458,7 +457,7 @@ protected: void constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, const MILU_VARIANT milu) const { - ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, milu, comm, amg ); + ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, milu, comm, amg ); } @@ -636,7 +635,7 @@ protected: { Vector weights(rhs_->size()); BlockVector rhs(0.0); - rhs[pressureIndex] = 1.0; + rhs[pressureVarIndex] = 1.0; int index = 0; ElementContext elemCtx(simulator_); const auto& vanguard = simulator_.vanguard(); @@ -657,7 +656,7 @@ protected: for (int ii = 0; ii < numEq; ++ii) { for (int jj = 0; jj < numEq; ++jj) { block[ii][jj] = storage[ii].derivative(jj)/storage_scale; - if (jj == 0) { + if (jj == pressureVarIndex) { block[ii][jj] *= pressure_scale; } } @@ -720,7 +719,7 @@ protected: Matrix& A = *matrix_; Vector weights(rhs_->size()); BlockVector rhs(0.0); - rhs[pressureIndex] = 1; + rhs[pressureVarIndex] = 1; const auto endi = A.end(); for (auto i = A.begin(); i!=endi; ++i) { const auto endj = (*i).end(); @@ -761,27 +760,20 @@ protected: const auto endj = (*i).end(); for (auto j = (*i).begin(); j != endj; ++j) { // assume it is something on all rows - // the blew logic depend on pressureIndex=0 Block& block = (*j); - for ( std::size_t ii = 0; ii < block.rows; ii++ ) { - if ( ii == 0 ) { - for (std::size_t jj = 0; jj < block.cols; jj++) { - block[0][jj] *= bweights[ii]; - } - } else { - for (std::size_t jj = 0; jj < block.cols; jj++) { - block[0][jj] += bweights[ii]*block[ii][jj]; - } + BlockVector neweq(0.0); + for (std::size_t ii = 0; ii < block.rows; ii++) { + for (std::size_t jj = 0; jj < block.cols; jj++) { + neweq[jj] += bweights[ii]*block[ii][jj]; } } + block[pressureEqnIndex] = neweq; } + BlockVector newrhs(0.0); for (std::size_t ii = 0; ii < brhs.size(); ii++) { - if ( ii == 0 ){ - brhs[0] *= bweights[ii]; - } else { - brhs[0] += bweights[ii]*brhs[ii]; - } + newrhs += bweights[ii]*brhs[ii]; } + brhs = newrhs; } } diff --git a/tests/test_blackoil_amg.cpp b/tests/test_blackoil_amg.cpp index 4fcbe7a9c..298425341 100644 --- a/tests/test_blackoil_amg.cpp +++ b/tests/test_blackoil_amg.cpp @@ -303,7 +303,7 @@ void runBlackoilAmgLaplace() smootherArgs.iterations = 1; Opm::CPRParameter param; - Opm::BlackoilAmg amg(param, + Opm::BlackoilAmg amg(param, {}, fop, criterion, smootherArgs, From b27e56e3864c5a480abdf067bf94743fc0540ee5 Mon Sep 17 00:00:00 2001 From: hnil Date: Mon, 18 Mar 2019 15:34:49 +0100 Subject: [PATCH 14/17] fixed bug in scaling of rhs for cpr --- opm/autodiff/ISTLSolverEbos.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index 3440956d3..b92991835 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -769,11 +769,11 @@ protected: } block[pressureEqnIndex] = neweq; } - BlockVector newrhs(0.0); + Scalar newrhs(0.0); for (std::size_t ii = 0; ii < brhs.size(); ii++) { newrhs += bweights[ii]*brhs[ii]; } - brhs = newrhs; + brhs[pressureEqnIndex] = newrhs; } } From 41acbcf861b837f22da91d2d99bd1c1a17b78ca4 Mon Sep 17 00:00:00 2001 From: hnil Date: Mon, 18 Mar 2019 16:11:04 +0100 Subject: [PATCH 15/17] change verbosity level setting for cpr to be more finegrained --- opm/autodiff/BlackoilAmg.hpp | 7 +++++-- opm/autodiff/CPRPreconditioner.hpp | 6 ++++-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index 88b3453a8..ecc9a2329 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -427,8 +427,11 @@ private: // Linear solver parameters const double tolerance = param_->cpr_solver_tol_; const int maxit = param_->cpr_max_iter_; - const int verbosity = ( param_->cpr_solver_verbose_ && - comm_.communicator().rank()==0 ) ? 1 : 0; + int verbosity = 0; + if (comm_.communicator().rank() == 0) { + verbosity = param_->cpr_solver_verbose_; + } + if ( param_->cpr_ell_solvetype_ == 0) { #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index 5ea90e94a..5a97aa10c 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -197,8 +197,10 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, const Vector& weights) { using AMG = BlackoilAmg; - const int verbosity = ( params.cpr_solver_verbose_ && comm.communicator().rank() == 0 ) ? 1 : 0; - + int verbosity = 0; + if (comm.communicator().rank() == 0) { + verbosity = params.cpr_solver_verbose_; + } // TODO: revise choice of parameters int coarsenTarget=1200; using Criterion = C; From 685a193558da9cc6fb61297cbd9f0a26d480b996 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 19 Mar 2019 15:06:55 +0100 Subject: [PATCH 16/17] make the branch self contained Concretely this avoids having to patch eWoms by adding a generic `Opm::transposeDenseMatrix()` template function instead of relying on the dense matrix class to provide a `transpose()` method. --- opm/autodiff/ISTLSolverEbos.hpp | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index b92991835..c1dccd9f3 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -61,6 +61,17 @@ END_PROPERTIES namespace Opm { +template +DenseMatrix transposeDenseMatrix(const DenseMatrix& M) +{ + DenseMatrix tmp; + for (int i = 0; i < M.rows; ++i) + for (int j = 0; j < M.cols; ++j) + tmp[j][i] = M[i][j]; + + return tmp; +} + //===================================================================== // Implementation for ISTL-matrix based operator //===================================================================== @@ -662,7 +673,7 @@ protected: } } BlockVector bweights; - MatrixBlockType block_transpose = block.transpose(); + MatrixBlockType block_transpose = Opm::transposeDenseMatrix(block); block_transpose.solve(bweights, rhs); bweights /= 1000.0; // given normal densities this scales weights to about 1. weights[index] = bweights; @@ -731,7 +742,7 @@ protected: } } BlockVector bweights; - auto diag_block_transpose = diag_block.transpose(); + auto diag_block_transpose = Opm::transposeDenseMatrix(diag_block); diag_block_transpose.solve(bweights, rhs); double abs_max = *std::max_element(bweights.begin(), bweights.end(), [](double a, double b){ return std::abs(a) < std::abs(b); } ); From 19b78b78c0fc967f4678465c754a7ef21b92baef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 20 Mar 2019 10:24:44 +0100 Subject: [PATCH 17/17] Updated after review comments. --- opm/autodiff/BlackoilAmg.hpp | 23 ++++++++++++++------- opm/autodiff/FlowLinearSolverParameters.hpp | 20 ++++++++---------- opm/autodiff/ISTLSolverEbos.hpp | 18 +++++++++++++--- 3 files changed, 39 insertions(+), 22 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index ecc9a2329..5f60c8528 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -426,7 +426,7 @@ private: } // Linear solver parameters const double tolerance = param_->cpr_solver_tol_; - const int maxit = param_->cpr_max_iter_; + const int maxit = param_->cpr_max_ell_iter_; int verbosity = 0; if (comm_.communicator().rank() == 0) { verbosity = param_->cpr_solver_verbose_; @@ -507,6 +507,12 @@ private: #endif } + // Warn if unknown options. + if (param_->cpr_ell_solvetype_ > 2 && comm_.communicator().rank() == 0) { + OpmLog::warning("cpr_ell_solver_type_unknown", "Unknown CPR elliptic solver type specification, using LoopSolver."); + } + + #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) delete sp; #endif @@ -926,16 +932,17 @@ private: * \brief An algebraic twolevel or multigrid approach for solving blackoil (supports CPR with and without AMG) * * This preconditioner first decouples the component used for coarsening using a simple scaling - * approach (e.g. Scheichl, Masson 2013,\see scaleMatrixDRS). Then it constructs the first - * coarse level system, either by simply extracting the coupling between the components at COMPONENT_INDEX - * in the matrix blocks or by extracting them and applying aggregation to them directly. This coarse level + * approach (e.g. Scheichl, Masson 2013,\see scaleMatrixDRS). Then it constructs the + * coarse level system. The coupling is defined by the weights corresponding to the element located at + * (COMPONENT_INDEX, VARIABLE_INDEX) in the block matrix. Then the coarse level system is constructed + * either by extracting these elements, or by applying aggregation to them directly. This coarse level * can be solved either by AMG or by ILU. The preconditioner is configured using CPRParameter. * \tparam O The type of the operator (encapsulating a BCRSMatrix). * \tparam S The type of the smoother. * \tparam C The type of coarsening criterion to use. * \tparam P The type of the class describing the parallelization. - * \tparam COMPONENT_INDEX The index of the component to use for coarsening (usually the pressure). - * \tparam VARIABLE_INDEX The index of the variable to use for coarsening (usually the pressure). + * \tparam COMPONENT_INDEX The index of the component to use for coarsening (usually water). + * \tparam VARIABLE_INDEX The index of the variable to use for coarsening (usually pressure). */ template @@ -1005,9 +1012,9 @@ public: : param_(param), 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_, diff --git a/opm/autodiff/FlowLinearSolverParameters.hpp b/opm/autodiff/FlowLinearSolverParameters.hpp index feeb77d7f..45b7a4339 100644 --- a/opm/autodiff/FlowLinearSolverParameters.hpp +++ b/opm/autodiff/FlowLinearSolverParameters.hpp @@ -63,7 +63,7 @@ NEW_PROP_TAG(SystemStrategy); NEW_PROP_TAG(ScaleLinearSystem); NEW_PROP_TAG(CprSolverVerbose); NEW_PROP_TAG(CprUseDrs); -NEW_PROP_TAG(CprMaxIter); +NEW_PROP_TAG(CprMaxEllIter); NEW_PROP_TAG(CprEllSolvetype); NEW_PROP_TAG(CprReuseSetup); @@ -83,11 +83,11 @@ SET_BOOL_PROP(FlowIstlSolverParams, UseAmg, false); SET_BOOL_PROP(FlowIstlSolverParams, UseCpr, false); SET_TYPE_PROP(FlowIstlSolverParams, LinearSolverBackend, Opm::ISTLSolverEbos); SET_BOOL_PROP(FlowIstlSolverParams, PreconditionerAddWellContributions, false); -SET_STRING_PROP(FlowIstlSolverParams, SystemStrategy, "original"); +SET_STRING_PROP(FlowIstlSolverParams, SystemStrategy, "none"); SET_BOOL_PROP(FlowIstlSolverParams, ScaleLinearSystem, false); SET_INT_PROP(FlowIstlSolverParams, CprSolverVerbose, 0); SET_BOOL_PROP(FlowIstlSolverParams, CprUseDrs, false); -SET_INT_PROP(FlowIstlSolverParams, CprMaxIter, 20); +SET_INT_PROP(FlowIstlSolverParams, CprMaxEllIter, 20); SET_INT_PROP(FlowIstlSolverParams, CprEllSolvetype, 0); SET_INT_PROP(FlowIstlSolverParams, CprReuseSetup, 0); @@ -111,7 +111,6 @@ namespace Opm bool cpr_ilu_reorder_sphere_; bool cpr_use_drs_; int cpr_max_ell_iter_; - int cpr_max_iter_; int cpr_ell_solvetype_; bool cpr_use_amg_; bool cpr_use_bicgstab_; @@ -130,7 +129,6 @@ namespace Opm cpr_max_ell_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; @@ -185,7 +183,7 @@ namespace Opm 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_max_ell_iter_ = EWOMS_GET_PARAM(TypeTag, int, CprMaxEllIter); cpr_ell_solvetype_ = EWOMS_GET_PARAM(TypeTag, int, CprEllSolvetype); cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup); } @@ -207,12 +205,12 @@ 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, std::string, SystemStrategy, "Strategy for reformulating and scaling linear system (none: no scaling -- should not be used with CPR, original: use weights that are equivalent to no scaling -- should not be used with CPR, simple: form pressure equation as simple sum of conservation equations, quasiimpes: form pressure equation based on diagonal block, trueimpes: form pressure equation based on linearization of accumulation term)"); 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, CprSolverVerbose, "Verbosity of cpr solver (0: silent, 1: print summary of inner linear solver, 2: print extensive information about inner linear solve, including setup information)"); + EWOMS_REGISTER_PARAM(TypeTag, bool, CprUseDrs, "Use dynamic row sum using weights"); + EWOMS_REGISTER_PARAM(TypeTag, int, CprMaxEllIter, "MaxIterations of the elliptic pressure part of the cpr solver"); + EWOMS_REGISTER_PARAM(TypeTag, int, CprEllSolvetype, "Solver type of elliptic pressure solve (0: bicgstab, 1: cg, 2: only amg preconditioner)"); EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse Amg Setup"); } diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index c1dccd9f3..a336d92eb 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -249,13 +249,16 @@ protected: bvec[pressureEqnIndex] = 1; weights_ = getSimpleWeights(bvec); } else { + if (parameters_.system_strategy_ != "none") { + OpmLog::warning("unknown_system_strategy", "Unknown linear solver system strategy: '" + parameters_.system_strategy_ + "', applying 'none' strategy."); + } form_cpr = false; } if (parameters_.scale_linear_system_) { // also scale weights this->scaleEquationsAndVariables(weights_); } - if (form_cpr && not(parameters_.cpr_use_drs_)) { + if (form_cpr && !(parameters_.cpr_use_drs_)) { scaleMatrixAndRhs(weights_); } if (weights_.size() == 0) { @@ -302,7 +305,7 @@ 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()) ); } @@ -642,6 +645,8 @@ protected: } // Weights to make approximate pressure equations. + // Calculated from the storage terms (only) of the + // conservation equations, ignoring all other terms. Vector getStorageWeights() const { Vector weights(rhs_->size()); @@ -682,6 +687,13 @@ protected: return weights; } + // Interaction between the CPR weights (the function argument 'weights') + // and the variable and equation weights from + // simulator_.model().primaryVarWeight() and + // simulator_.model().eqWeight() is nontrivial and does not work + // at the moment. Possibly refactoring of ewoms weight treatment + // is needed. In the meantime this function shows what needs to be + // done to integrate the weights properly. void scaleEquationsAndVariables(Vector& weights) { // loop over primary variables @@ -702,7 +714,7 @@ protected: for (std::size_t ii = 0; ii < brhs.size(); ii++) { brhs[ii] *= simulator_.model().eqWeight(i.index(), ii); } - if (weights_.size() == matrix_->N()) { + if (weights.size() == matrix_->N()) { BlockVector& bw = weights[i.index()]; for (std::size_t ii = 0; ii < brhs.size(); ii++) { bw[ii] /= simulator_.model().eqWeight(i.index(), ii);