154 using Matrix =
typename SparseMatrixAdapter::IstlMatrix;
157 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
158 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
165 constexpr static bool isIncompatibleWithCprw = 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") {
223 "The polymer injectivity model is incompatible with the CPRW linear solver.\n"
224 "Choose a different option, for example --linear-solver=ilu0");
228 if (parameters_[0].linsolver_ ==
"hybrid") {
236 FlowLinearSolverParameters
para;
238 para.linsolver_ =
"cprw";
239 parameters_.push_back(
para);
241 Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
242 Parameters::IsSet<Parameters::LinearSolverReduction>()));
245 FlowLinearSolverParameters
para;
247 para.linsolver_ =
"ilu0";
248 parameters_.push_back(
para);
250 Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
251 Parameters::IsSet<Parameters::LinearSolverReduction>()));
255 assert(parameters_.size() == 1);
259 if (parameters_[0].is_nldd_local_solver_) {
261 Parameters::IsSet<Parameters::NlddLocalLinearSolverMaxIter>(),
262 Parameters::IsSet<Parameters::NlddLocalLinearSolverReduction>()));
266 Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
267 Parameters::IsSet<Parameters::LinearSolverReduction>()));
270 flexibleSolver_.resize(prm_.size());
272 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
274 comm_.reset(
new CommunicationType( simulator_.vanguard().grid().comm() ) );
281 ElementMapper
elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
282 detail::findOverlapAndInterior(simulator_.vanguard().grid(),
elemMapper, overlapRows_, interiorRows_);
283 useWellConn_ = Parameters::Get<Parameters::MatrixAddWellContributions>();
284 const bool ownersFirst = Parameters::Get<Parameters::OwnerCellsFirst>();
286 const std::string
msg =
"The linear solver no longer supports --owner-cells-first=false.";
293 const int interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(),
true);
294 for (
auto&
f : flexibleSolver_) {
295 f.interiorCellNum_ = interiorCellNum_;
300 const std::size_t size = simulator_.vanguard().grid().leafGridView().size(0);
301 detail::copyParValues(parallelInformation_, size, *comm_);
306 if (
on_io_rank && parameters_[activeSolverNum_].linear_solver_print_json_definition_) {
307 std::ostringstream
os;
308 os <<
"\nProperty tree for linear solvers:\n";
309 for (std::size_t i = 0; i<prm_.size(); i++) {
310 prm_[i].write_json(
os,
true);
312 OpmLog::note(
os.str());
321 void setActiveSolver(
const int num)
323 if (
num >
static_cast<int>(prm_.size()) - 1) {
324 OPM_THROW(std::logic_error,
"Solver number " + std::to_string(
num) +
" not available.");
326 activeSolverNum_ =
num;
327 if (simulator_.gridView().comm().rank() == 0) {
328 OpmLog::debug(
"Active solver = " + std::to_string(activeSolverNum_)
329 +
" (" + parameters_[activeSolverNum_].linsolver_ +
")");
333 int numAvailableSolvers()
335 return flexibleSolver_.size();
338 void initPrepare(
const Matrix& M, Vector&
b)
340 const bool firstcall = (matrix_ ==
nullptr);
347 matrix_ =
const_cast<Matrix*
>(&M);
349 useWellConn_ = Parameters::Get<Parameters::MatrixAddWellContributions>();
353 if ( &M != matrix_ ) {
355 "Matrix objects are expected to be reused when reassembling!");
361 if (isParallel() && prm_[activeSolverNum_].
template get<std::string>(
"preconditioner.type") !=
"ParOverILU0") {
362 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
366 void prepare(
const SparseMatrixAdapter& M, Vector&
b)
368 prepare(M.istlMatrix(),
b);
371 void prepare(
const Matrix& M, Vector&
b)
377 prepareFlexibleSolver();
382 void setResidual(Vector& )
387 void getResidual(Vector&
b)
const
392 void setMatrix(
const SparseMatrixAdapter& )
397 int getSolveCount()
const {
401 void resetSolveCount() {
405 bool solve(Vector& x)
410 const int verbosity = prm_[activeSolverNum_].get(
"verbosity", 0);
413 Helper::writeSystem(simulator_,
420 Dune::InverseOperatorResult
result;
423 assert(flexibleSolver_[activeSolverNum_].solver_);
424 flexibleSolver_[activeSolverNum_].solver_->apply(x, *rhs_,
result);
446 const CommunicationType* comm()
const {
return comm_.get(); }
448 void setDomainIndex(
const int index)
450 domainIndex_ = index;
453 bool isNlddLocalSolver()
const
455 return parameters_[activeSolverNum_].is_nldd_local_solver_;
460 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
463 void checkConvergence(
const Dune::InverseOperatorResult&
result )
const
466 iterations_ =
result.iterations;
467 converged_ =
result.converged;
469 if(
result.reduction < parameters_[activeSolverNum_].relaxed_linear_solver_reduction_){
470 std::stringstream
ss;
471 ss<<
"Full linear solver tolerance not achieved. The reduction is:" <<
result.reduction
472 <<
" after " <<
result.iterations <<
" iterations ";
473 OpmLog::warning(
ss.str());
478 if (!parameters_[activeSolverNum_].ignoreConvergenceFailure_ && !converged_) {
479 const std::string
msg(
"Convergence failure for linear solver.");
485 bool isParallel()
const {
487 return !forceSerial_ && comm_->communicator().size() > 1;
493 void prepareFlexibleSolver()
498 if (isNlddLocalSolver()) {
499 auto wellOp = std::make_unique<DomainWellModelAsLinearOperator<WellModel, Vector, Vector>>(simulator_.problem().wellModel());
500 wellOp->setDomainIndex(domainIndex_);
501 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(
wellOp);
504 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
505 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(
wellOp);
508 std::function<Vector()>
weightCalculator = this->getWeightsCalculator(prm_[activeSolverNum_], getMatrix(), pressureIndex);
510 flexibleSolver_[activeSolverNum_].create(getMatrix(),
512 prm_[activeSolverNum_],
521 flexibleSolver_[activeSolverNum_].pre_->update();
532 if (flexibleSolver_.empty()) {
535 if (!flexibleSolver_[activeSolverNum_].solver_) {
539 if (flexibleSolver_[activeSolverNum_].pre_->hasPerfectUpdate()) {
545 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 0) {
549 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 1) {
551 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
554 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 2) {
558 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 3) {
562 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 4) {
564 const int step = this->parameters_[activeSolverNum_].cpr_reuse_interval_;
565 const bool create = ((solveCount_ % step) == 0);
569 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
570 std::string
msg =
"Invalid value: " + std::to_string(this->parameters_[activeSolverNum_].cpr_reuse_setup_)
571 +
" for --cpr-reuse-setup parameter, run with --help to see allowed values.";
575 throw std::runtime_error(
msg);
584 std::function<Vector()> getWeightsCalculator(
const PropertyTree& prm,
585 const Matrix& matrix,
590 using namespace std::string_literals;
596 const auto weightsType = prm.
get(
"preconditioner.weight_type"s,
"quasiimpes"s);
602 return Amg::getQuasiImpesWeights<Matrix, Vector>(matrix,
610 Vector weights(rhs_->size());
611 ElementContext elemCtx(simulator_);
613 simulator_.vanguard().gridView(),
614 elemCtx, simulator_.model(),
622 Vector weights(rhs_->size());
623 ElementContext elemCtx(simulator_);
624 Amg::getTrueImpesWeightsAnalytic(
pressIndex, weights,
625 simulator_.vanguard().gridView(),
626 elemCtx, simulator_.model(),
633 "not implemented for cpr."
634 " Please use quasiimpes, trueimpes or trueimpesanalytic.");
646 const Matrix& getMatrix()
const
651 const Simulator& simulator_;
652 mutable int iterations_;
653 mutable int solveCount_;
654 mutable bool converged_;
655 std::any parallelInformation_;
661 int activeSolverNum_ = 0;
662 std::vector<detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType>> flexibleSolver_;
663 std::vector<int> overlapRows_;
664 std::vector<int> interiorRows_;
666 int domainIndex_ = -1;
670 std::vector<FlowLinearSolverParameters> parameters_;
671 bool forceSerial_ =
false;
672 std::vector<PropertyTree> prm_;
674 std::shared_ptr< CommunicationType > comm_;