153 using Matrix =
typename SparseMatrixAdapter::IstlMatrix;
156 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
157 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
165 constexpr static bool isIncompatibleWithCprw = enableMICP || enablePolymerMolarWeight;
168 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
170 using CommunicationType = Dune::Communication<int>;
174 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
176 static void registerParameters()
178 FlowLinearSolverParameters::registerParameters();
191 : simulator_(simulator),
204 : simulator_(simulator),
210 parameters_.resize(1);
211 parameters_[0].init(simulator_.vanguard().eclState().getSimulationConfig().useCPR());
219 if (isIncompatibleWithCprw) {
221 if (parameters_[0].linsolver_ ==
"cprw" || parameters_[0].linsolver_ ==
"hybrid") {
225 }
else if (enablePolymerMolarWeight) {
230 "Choose a different option, for example --linear-solver=ilu0");
234 if (parameters_[0].linsolver_ ==
"hybrid") {
242 FlowLinearSolverParameters
para;
244 para.linsolver_ =
"cprw";
245 parameters_.push_back(
para);
247 Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
248 Parameters::IsSet<Parameters::LinearSolverReduction>()));
251 FlowLinearSolverParameters
para;
253 para.linsolver_ =
"ilu0";
254 parameters_.push_back(
para);
256 Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
257 Parameters::IsSet<Parameters::LinearSolverReduction>()));
261 assert(parameters_.size() == 1);
265 if (parameters_[0].is_nldd_local_solver_) {
267 Parameters::IsSet<Parameters::NlddLocalLinearSolverMaxIter>(),
268 Parameters::IsSet<Parameters::NlddLocalLinearSolverReduction>()));
272 Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
273 Parameters::IsSet<Parameters::LinearSolverReduction>()));
276 flexibleSolver_.resize(prm_.size());
278 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
280 comm_.reset(
new CommunicationType( simulator_.vanguard().grid().comm() ) );
287 ElementMapper
elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
288 detail::findOverlapAndInterior(simulator_.vanguard().grid(),
elemMapper, overlapRows_, interiorRows_);
289 useWellConn_ = Parameters::Get<Parameters::MatrixAddWellContributions>();
290 const bool ownersFirst = Parameters::Get<Parameters::OwnerCellsFirst>();
292 const std::string
msg =
"The linear solver no longer supports --owner-cells-first=false.";
299 const int interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(),
true);
300 for (
auto&
f : flexibleSolver_) {
301 f.interiorCellNum_ = interiorCellNum_;
306 const std::size_t size = simulator_.vanguard().grid().leafGridView().size(0);
307 detail::copyParValues(parallelInformation_, size, *comm_);
312 if (
on_io_rank && parameters_[activeSolverNum_].linear_solver_print_json_definition_) {
313 std::ostringstream
os;
314 os <<
"Property tree for linear solvers:\n";
315 for (std::size_t i = 0; i<prm_.size(); i++) {
316 prm_[i].write_json(
os,
true);
318 OpmLog::note(
os.str());
327 void setActiveSolver(
const int num)
329 if (
num >
static_cast<int>(prm_.size()) - 1) {
330 OPM_THROW(std::logic_error,
"Solver number " + std::to_string(
num) +
" not available.");
332 activeSolverNum_ =
num;
333 if (simulator_.gridView().comm().rank() == 0) {
334 OpmLog::debug(
"Active solver = " + std::to_string(activeSolverNum_)
335 +
" (" + parameters_[activeSolverNum_].linsolver_ +
")");
339 int numAvailableSolvers()
341 return flexibleSolver_.size();
344 void initPrepare(
const Matrix& M, Vector&
b)
346 const bool firstcall = (matrix_ ==
nullptr);
353 matrix_ =
const_cast<Matrix*
>(&M);
355 useWellConn_ = Parameters::Get<Parameters::MatrixAddWellContributions>();
359 if ( &M != matrix_ ) {
361 "Matrix objects are expected to be reused when reassembling!");
367 if (isParallel() && prm_[activeSolverNum_].
template get<std::string>(
"preconditioner.type") !=
"ParOverILU0") {
368 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
372 void prepare(
const SparseMatrixAdapter& M, Vector&
b)
374 prepare(M.istlMatrix(),
b);
377 void prepare(
const Matrix& M, Vector&
b)
383 prepareFlexibleSolver();
387 void setResidual(Vector& )
392 void getResidual(Vector&
b)
const
397 void setMatrix(
const SparseMatrixAdapter& )
402 int getSolveCount()
const {
406 void resetSolveCount() {
410 bool solve(Vector& x)
415 const int verbosity = prm_[activeSolverNum_].get(
"verbosity", 0);
418 Helper::writeSystem(simulator_,
425 Dune::InverseOperatorResult
result;
428 assert(flexibleSolver_[activeSolverNum_].solver_);
429 flexibleSolver_[activeSolverNum_].solver_->apply(x, *rhs_,
result);
451 const CommunicationType* comm()
const {
return comm_.get(); }
453 void setDomainIndex(
const int index)
455 domainIndex_ = index;
458 bool isNlddLocalSolver()
const
460 return parameters_[activeSolverNum_].is_nldd_local_solver_;
465 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
468 void checkConvergence(
const Dune::InverseOperatorResult&
result )
const
471 iterations_ =
result.iterations;
472 converged_ =
result.converged;
474 if(
result.reduction < parameters_[activeSolverNum_].relaxed_linear_solver_reduction_){
475 std::stringstream
ss;
476 ss<<
"Full linear solver tolerance not achieved. The reduction is:" <<
result.reduction
477 <<
" after " <<
result.iterations <<
" iterations ";
478 OpmLog::warning(
ss.str());
483 if (!parameters_[activeSolverNum_].ignoreConvergenceFailure_ && !converged_) {
484 const std::string
msg(
"Convergence failure for linear solver.");
490 bool isParallel()
const {
492 return !forceSerial_ && comm_->communicator().size() > 1;
498 void prepareFlexibleSolver()
503 if (isNlddLocalSolver()) {
504 auto wellOp = std::make_unique<DomainWellModelAsLinearOperator<WellModel, Vector, Vector>>(simulator_.problem().wellModel());
505 wellOp->setDomainIndex(domainIndex_);
506 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(
wellOp);
509 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
510 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(
wellOp);
513 std::function<Vector()>
weightCalculator = this->getWeightsCalculator(prm_[activeSolverNum_], getMatrix(), pressureIndex);
515 flexibleSolver_[activeSolverNum_].create(getMatrix(),
517 prm_[activeSolverNum_],
526 flexibleSolver_[activeSolverNum_].pre_->update();
537 if (flexibleSolver_.empty()) {
540 if (!flexibleSolver_[activeSolverNum_].solver_) {
544 if (flexibleSolver_[activeSolverNum_].pre_->hasPerfectUpdate()) {
550 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 0) {
554 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 1) {
556 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
559 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 2) {
563 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 3) {
567 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 4) {
569 const int step = this->parameters_[activeSolverNum_].cpr_reuse_interval_;
570 const bool create = ((solveCount_ % step) == 0);
574 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
575 std::string
msg =
"Invalid value: " + std::to_string(this->parameters_[activeSolverNum_].cpr_reuse_setup_)
576 +
" for --cpr-reuse-setup parameter, run with --help to see allowed values.";
580 throw std::runtime_error(
msg);
589 std::function<Vector()> getWeightsCalculator(
const PropertyTree& prm,
590 const Matrix& matrix,
595 using namespace std::string_literals;
601 const auto weightsType = prm.get(
"preconditioner.weight_type"s,
"quasiimpes"s);
607 return Amg::getQuasiImpesWeights<Matrix, Vector>(matrix,
615 Vector weights(rhs_->size());
616 ElementContext
elemCtx(simulator_);
618 simulator_.vanguard().gridView(),
627 Vector weights(rhs_->size());
628 ElementContext
elemCtx(simulator_);
629 Amg::getTrueImpesWeightsAnalytic(
pressIndex, weights,
630 simulator_.vanguard().gridView(),
638 "not implemented for cpr."
639 " Please use quasiimpes, trueimpes or trueimpesanalytic.");
651 const Matrix& getMatrix()
const
656 const Simulator& simulator_;
657 mutable int iterations_;
658 mutable int solveCount_;
659 mutable bool converged_;
660 std::any parallelInformation_;
666 int activeSolverNum_ = 0;
667 std::vector<detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType>> flexibleSolver_;
668 std::vector<int> overlapRows_;
669 std::vector<int> interiorRows_;
671 int domainIndex_ = -1;
675 std::vector<FlowLinearSolverParameters> parameters_;
676 bool forceSerial_ =
false;
677 std::vector<PropertyTree> prm_;
679 std::shared_ptr< CommunicationType > comm_;