Merge pull request #6 from akva2/improve_cpr_mod

Improve cpr mod
This commit is contained in:
Atgeirr Flø Rasmussen 2019-03-14 11:53:16 +01:00 committed by GitHub
commit 81ba405f26
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 53 additions and 66 deletions

View File

@ -79,7 +79,6 @@ std::unique_ptr< Dune::MatrixAdapter<M,X,Y> > createOperatorPtr(const Dune::Matr
{ {
return std::make_unique< Dune::MatrixAdapter<M,X,Y> >(matrix); return std::make_unique< Dune::MatrixAdapter<M,X,Y> >(matrix);
} }
/** /**
* \brief Creates an OverlappingSchwarzOperator as an operator. * \brief Creates an OverlappingSchwarzOperator as an operator.
* *
@ -119,13 +118,9 @@ scaleMatrixDRS(const Operator& op, const Communication& comm,
for (auto j = (*i).begin(); j != endj; ++j) { for (auto j = (*i).begin(); j != endj; ++j) {
BlockVector bvec(0.0); BlockVector bvec(0.0);
Block& block = *j; Block& block = *j;
for (std::size_t ii = 0; ii < Block::rows; ii++) { BlockVector& bvec = block[pressureIndex];
for (std::size_t jj = 0; jj < Block::cols; jj++) { // should introduce limits which also change the weights
// should introduce limmits which also change the weights block.mtv(bw, bvec);
bvec[jj] += bw[ii]*block[ii][jj];
}
}
block[pressureIndex] = bvec;
} }
} }
} }
@ -144,13 +139,9 @@ void scaleVectorDRS(Vector& vector, std::size_t pressureIndex, const Opm::CPRPar
using Block = typename Vector::block_type; using Block = typename Vector::block_type;
if (param.cpr_use_drs_) { if (param.cpr_use_drs_) {
for (std::size_t j = 0; j < vector.size(); ++j) { for (std::size_t j = 0; j < vector.size(); ++j) {
double val(0.0);
Block& block = vector[j]; Block& block = vector[j];
const Block& bw = weights[j]; const Block& bw = weights[j];
for (std::size_t i = 0; i < Block::dimension; i++) { block[pressureIndex] = bw.dot(block);
val += bw[i]*block[i];
}
block[pressureIndex] = val;
} }
} }
} }
@ -995,8 +986,7 @@ public:
levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_), levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_),
coarseSolverPolicy_(&param, smargs, criterion), coarseSolverPolicy_(&param, smargs, criterion),
twoLevelMethod_(std::get<1>(scaledMatrixOperator_), smoother_, twoLevelMethod_(std::get<1>(scaledMatrixOperator_), smoother_,
levelTransferPolicy_, levelTransferPolicy_, coarseSolverPolicy_, 0, 1)
coarseSolverPolicy_, 0, 1)
{} {}
void pre(typename TwoLevelMethod::FineDomainType& x, void pre(typename TwoLevelMethod::FineDomainType& x,

View File

@ -188,7 +188,6 @@ namespace Opm
cpr_max_iter_ = EWOMS_GET_PARAM(TypeTag, int, CprMaxIter); cpr_max_iter_ = EWOMS_GET_PARAM(TypeTag, int, CprMaxIter);
cpr_ell_solvetype_ = EWOMS_GET_PARAM(TypeTag, int, CprEllSolvetype); cpr_ell_solvetype_ = EWOMS_GET_PARAM(TypeTag, int, CprEllSolvetype);
cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup); cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup);
} }
template <class TypeTag> template <class TypeTag>
@ -215,7 +214,6 @@ namespace Opm
EWOMS_REGISTER_PARAM(TypeTag, int, CprMaxIter, "MaxIterations of the pressure amg solver"); 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, 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, int, CprReuseSetup, "Reuse Amg Setup");
} }
FlowLinearSolverParameters() { reset(); } FlowLinearSolverParameters() { reset(); }

View File

@ -253,6 +253,10 @@ protected:
parameters_.cpr_use_drs_ = false; parameters_.cpr_use_drs_ = false;
} }
} else { } 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_) { if (parameters_.scale_linear_system_) {
// also scale weights // also scale weights
this->scaleEquationsAndVariables(weights_); this->scaleEquationsAndVariables(weights_);
@ -294,6 +298,7 @@ protected:
} }
else else
{ {
const WellModel& wellModel = simulator_.problem().wellModel();
typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false > Operator; typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false > Operator;
Operator opA(*matrix_, *matrix_, wellModel); Operator opA(*matrix_, *matrix_, wellModel);
solve( opA, x, *rhs_ ); solve( opA, x, *rhs_ );
@ -304,7 +309,6 @@ protected:
} }
return converged_; return converged_;
} }

View File

@ -303,13 +303,8 @@ void runBlackoilAmgLaplace()
smootherArgs.iterations = 1; smootherArgs.iterations = 1;
Opm::CPRParameter param; Opm::CPRParameter param;
Vector weights(b.size());
for (auto& elem : weights) {
elem = 1.0;
}
Opm::BlackoilAmg<Operator,ParSmoother,Criterion,Communication,0> amg(param, Opm::BlackoilAmg<Operator,ParSmoother,Criterion,Communication,0> amg(param,
weights, {},
fop, criterion, fop, criterion,
smootherArgs, smootherArgs,
comm); comm);