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

Changes to make cpr work

Browse Source
This commit is contained in:
hnil 2019-03-12 13:55:11 +01:00
parent 6452fbc3cd
commit 40537f1999
4 changed files with 491 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

114
opm/autodiff/BlackoilAmg.hpp
View File

@ -21,6 +21,7 @@
#include <ewoms/linear/matrixblock.hh>
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
#include <opm/autodiff/FlowLinearSolverParameters.hpp>
#include <opm/autodiff/CPRPreconditioner.hpp>
#include <dune/istl/paamg/twolevelmethod.hh>
#include <dune/istl/paamg/aggregates.hh>
@ -67,6 +68,11 @@ Dune::MatrixAdapter<M,X,Y> createOperator(const Dune::MatrixAdapter<M,X,Y>&, con
return Dune::MatrixAdapter<M,X,Y>(matrix);
}
template<class M, class X, class Y, class T>
std::unique_ptr< Dune::MatrixAdapter<M,X,Y> > createOperatorPtr(const Dune::MatrixAdapter<M,X,Y>&, const M& matrix, const T&)
{
return std::make_unique< Dune::MatrixAdapter<M,X,Y> >(matrix);
}
/**
* \brief Creates an OverlappingSchwarzOperator as an operator.
*
@ -89,30 +95,35 @@ Dune::OverlappingSchwarzOperator<M,X,Y,T> 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<class Operator, class Communication>
template<class Operator, class Communication,class Vector>
std::tuple<std::unique_ptr<typename Operator::matrix_type>, 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> 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<class Vector>
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<X> solver(const_cast<typename AMGType::Operator&>(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<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp, *prec,
@ -453,6 +465,30 @@ private:
}
#endif
}
else
{
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
Dune::LoopSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp, *prec,
tolerance, maxit, verbosity);
solver.apply(x,b,res);
#else
if ( !amg_ )
{
Dune::LoopSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp,
reinterpret_cast<Smoother&>(*prec),
tolerance, maxit, verbosity);
solver.apply(x,b,res);
}
else
{
Dune::LoopSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp,
reinterpret_cast<AMGType&>(*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<Smoother>(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<std::unique_ptr<Matrix>, Operator> scaledMatrixOperator_;
std::shared_ptr<Smoother> smoother_;
LevelTransferPolicy levelTransferPolicy_;
CoarseSolverPolicy coarseSolverPolicy_;
TwoLevelMethod twoLevelMethod_;
//BlockVector weights_;
};
namespace ISTLUtility

12
opm/autodiff/CPRPreconditioner.hpp
View File

@ -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<Op,S,C,P,index> >& amgPtr,
const CPRParameter& params)
const CPRParameter& params,
const Vector& weights)
{
using AMG = BlackoilAmg<Op,S,C,P,index>;
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 >

136
opm/autodiff/FlowLinearSolverParameters.hpp
View File

@ -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<TypeTag>);
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 <class TypeTag>
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 <class TypeTag>
@ -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()

331
opm/autodiff/ISTLSolverEbos.hpp
View File

@ -34,6 +34,7 @@
#include <ewoms/common/parametersystem.hh>
#include <ewoms/common/propertysystem.hh>
//#include <ewoms/linear/matrixmarket_ewoms.hh>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/operators.hh>
@ -170,6 +171,7 @@ protected:
template <class TypeTag>
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<C>( *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</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
Dune::FieldVector<Evaluation, numEq> 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> matrix_;
Vector *rhs_;
std::unique_ptr<Matrix> matrix_for_preconditioner_;
std::vector<std::pair<int,std::vector<int>>> overlapRowAndColumns_;
FlowLinearSolverParameters parameters_;
Vector weights_;
bool scale_variables_;
}; // end ISTLSolver
} // namespace Opm