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();
190 bool forceSerial =
false)
191 : simulator_(simulator),
195 parameters_{parameters},
196 forceSerial_(forceSerial)
204 : simulator_(simulator),
210 parameters_.resize(1);
211 parameters_[0].init(simulator_.vanguard().eclState().getSimulationConfig().useCPR());
217 OPM_TIMEBLOCK(IstlSolver);
219 if (isIncompatibleWithCprw) {
221 if (parameters_[0].linsolver_ ==
"cprw" || parameters_[0].linsolver_ ==
"hybrid") {
222 std::string incompatible_model =
"Unknown";
224 incompatible_model =
"MICP";
225 }
else if (enablePolymerMolarWeight) {
226 incompatible_model =
"Polymer injectivity";
228 OPM_THROW(std::runtime_error,
229 incompatible_model +
" model is incompatible with the CPRW linear solver.\n"
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);
262 assert(prm_.empty());
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() ) );
282 extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
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.";
296 OPM_THROW_NOLOG(std::runtime_error, msg);
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_ ) {
360 OPM_THROW(std::logic_error,
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)
379 OPM_TIMEBLOCK(istlSolverPrepare);
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)
412 OPM_TIMEBLOCK(istlSolverSolve);
415 const int verbosity = prm_[activeSolverNum_].get(
"verbosity", 0);
416 const bool write_matrix = verbosity > 10;
418 Helper::writeSystem(simulator_,
425 Dune::InverseOperatorResult result;
427 OPM_TIMEBLOCK(flexibleSolverApply);
428 assert(flexibleSolver_[activeSolverNum_].solver_);
429 flexibleSolver_[activeSolverNum_].solver_->apply(x, *rhs_, result);
433 checkConvergence(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.");
485 OPM_THROW_NOLOG(NumericalProblem, msg);
490 bool isParallel()
const {
492 return !forceSerial_ && comm_->communicator().size() > 1;
498 void prepareFlexibleSolver()
500 OPM_TIMEBLOCK(flexibleSolverPrepare);
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);
514 OPM_TIMEBLOCK(flexibleSolverCreate);
515 flexibleSolver_[activeSolverNum_].create(getMatrix(),
517 prm_[activeSolverNum_],
525 OPM_TIMEBLOCK(flexibleSolverUpdate);
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();
557 return newton_iteration == 0;
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,
591 std::size_t pressIndex)
const
593 std::function<Vector()> weightsCalculator;
595 using namespace std::string_literals;
597 auto preconditionerType = prm.get(
"preconditioner.type"s,
"cpr"s);
598 if (preconditionerType ==
"cpr" || preconditionerType ==
"cprt"
599 || preconditionerType ==
"cprw" || preconditionerType ==
"cprwt") {
600 const bool transpose = preconditionerType ==
"cprt" || preconditionerType ==
"cprwt";
601 const auto weightsType = prm.get(
"preconditioner.weight_type"s,
"quasiimpes"s);
602 if (weightsType ==
"quasiimpes") {
606 weightsCalculator = [matrix, transpose, pressIndex]() {
607 return Amg::getQuasiImpesWeights<Matrix, Vector>(matrix,
611 }
else if ( weightsType ==
"trueimpes" ) {
615 Vector weights(rhs_->size());
616 ElementContext elemCtx(simulator_);
617 Amg::getTrueImpesWeights(pressIndex, weights,
618 simulator_.vanguard().gridView(),
619 elemCtx, simulator_.model(),
623 }
else if (weightsType ==
"trueimpesanalytic" ) {
627 Vector weights(rhs_->size());
628 ElementContext elemCtx(simulator_);
629 Amg::getTrueImpesWeightsAnalytic(pressIndex, weights,
630 simulator_.vanguard().gridView(),
631 elemCtx, simulator_.model(),
636 OPM_THROW(std::invalid_argument,
637 "Weights type " + weightsType +
638 "not implemented for cpr."
639 " Please use quasiimpes, trueimpes or trueimpesanalytic.");
642 return weightsCalculator;
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_;