From 40537f199914715ee5b03829c09ec213fca3c996 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 12 Mar 2019 13:55:11 +0100 Subject: [PATCH] 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