145 using GridView = GetPropType<TypeTag, Properties::GridView>;
146 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
147 using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
148 using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
149 using Indices = GetPropType<TypeTag, Properties::Indices>;
150 using WellModel = GetPropType<TypeTag, Properties::WellModel>;
151 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
152 using Matrix =
typename SparseMatrixAdapter::IstlMatrix;
153 using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
154 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
155 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
156 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
159 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
160 constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
163 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
165 using CommunicationType = Dune::CollectiveCommunication<int>;
169 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
171 static void registerParameters()
173 FlowLinearSolverParameters::registerParameters<TypeTag>();
185 bool forceSerial =
false)
186 : simulator_(simulator),
190 parameters_{parameters},
191 forceSerial_(forceSerial)
199 : simulator_(simulator),
205 parameters_.resize(1);
206 parameters_[0].template init<TypeTag>(simulator_.vanguard().eclState().getSimulationConfig().useCPR());
212 OPM_TIMEBLOCK(IstlSolver);
214 if (parameters_[0].linsolver_ ==
"hybrid") {
223 para.init<TypeTag>(
false);
224 para.linsolver_ =
"cprw";
225 parameters_.push_back(para);
227 Parameters::isSet<TypeTag,Properties::LinearSolverMaxIter>(),
228 Parameters::isSet<TypeTag,Properties::LinearSolverReduction>()));
231 FlowLinearSolverParameters para;
232 para.init<TypeTag>(
false);
233 para.linsolver_ =
"ilu0";
234 parameters_.push_back(para);
236 Parameters::isSet<TypeTag,Properties::LinearSolverMaxIter>(),
237 Parameters::isSet<TypeTag,Properties::LinearSolverReduction>()));
242 assert(parameters_.size() == 1);
243 assert(prm_.empty());
245 Parameters::isSet<TypeTag,Properties::LinearSolverMaxIter>(),
246 Parameters::isSet<TypeTag,Properties::LinearSolverReduction>()));
248 flexibleSolver_.resize(prm_.size());
250 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
252 comm_.reset(
new CommunicationType( simulator_.vanguard().grid().comm() ) );
254 extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
259 ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
260 detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
261 useWellConn_ = Parameters::get<TypeTag, Properties::MatrixAddWellContributions>();
262 const bool ownersFirst = Parameters::get<TypeTag, Properties::OwnerCellsFirst>();
264 const std::string msg =
"The linear solver no longer supports --owner-cells-first=false.";
268 OPM_THROW_NOLOG(std::runtime_error, msg);
271 const int interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(),
true);
272 for (
auto& f : flexibleSolver_) {
273 f.interiorCellNum_ = interiorCellNum_;
278 const std::size_t size = simulator_.vanguard().grid().leafGridView().size(0);
279 detail::copyParValues(parallelInformation_, size, *comm_);
284 if (on_io_rank && parameters_[activeSolverNum_].linear_solver_print_json_definition_) {
285 std::ostringstream os;
286 os <<
"Property tree for linear solvers:\n";
287 for (std::size_t i = 0; i<prm_.size(); i++) {
288 prm_[i].write_json(os,
true);
290 OpmLog::note(os.str());
299 void setActiveSolver(
const int num)
301 if (num >
static_cast<int>(prm_.size()) - 1) {
302 OPM_THROW(std::logic_error,
"Solver number " + std::to_string(num) +
" not available.");
304 activeSolverNum_ = num;
305 if (simulator_.gridView().comm().rank() == 0) {
306 OpmLog::debug(
"Active solver = " + std::to_string(activeSolverNum_)
307 +
" (" + parameters_[activeSolverNum_].linsolver_ +
")");
311 int numAvailableSolvers()
313 return flexibleSolver_.size();
316 void initPrepare(
const Matrix& M, Vector& b)
318 const bool firstcall = (matrix_ ==
nullptr);
325 matrix_ =
const_cast<Matrix*
>(&M);
327 useWellConn_ = Parameters::get<TypeTag, Properties::MatrixAddWellContributions>();
331 if ( &M != matrix_ ) {
332 OPM_THROW(std::logic_error,
333 "Matrix objects are expected to be reused when reassembling!");
339 if (isParallel() && prm_[activeSolverNum_].
template get<std::string>(
"preconditioner.type") !=
"ParOverILU0") {
340 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
344 void prepare(
const SparseMatrixAdapter& M, Vector& b)
346 prepare(M.istlMatrix(), b);
349 void prepare(
const Matrix& M, Vector& b)
351 OPM_TIMEBLOCK(istlSolverPrepare);
355 prepareFlexibleSolver();
359 void setResidual(Vector& )
364 void getResidual(Vector& b)
const
369 void setMatrix(
const SparseMatrixAdapter& )
374 int getSolveCount()
const {
378 void resetSolveCount() {
382 bool solve(Vector& x)
384 OPM_TIMEBLOCK(istlSolverSolve);
387 const int verbosity = prm_[activeSolverNum_].get(
"verbosity", 0);
388 const bool write_matrix = verbosity > 10;
390 Helper::writeSystem(simulator_,
397 Dune::InverseOperatorResult result;
399 OPM_TIMEBLOCK(flexibleSolverApply);
400 assert(flexibleSolver_[activeSolverNum_].solver_);
401 flexibleSolver_[activeSolverNum_].solver_->apply(x, *rhs_, result);
405 checkConvergence(result);
423 const CommunicationType* comm()
const {
return comm_.get(); }
427 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
430 void checkConvergence(
const Dune::InverseOperatorResult& result )
const
433 iterations_ = result.iterations;
434 converged_ = result.converged;
436 if(result.reduction < parameters_[activeSolverNum_].relaxed_linear_solver_reduction_){
437 std::stringstream ss;
438 ss<<
"Full linear solver tolerance not achieved. The reduction is:" << result.reduction
439 <<
" after " << result.iterations <<
" iterations ";
440 OpmLog::warning(ss.str());
445 if (!parameters_[activeSolverNum_].ignoreConvergenceFailure_ && !converged_) {
446 const std::string msg(
"Convergence failure for linear solver.");
447 OPM_THROW_NOLOG(NumericalProblem, msg);
452 bool isParallel()
const {
454 return !forceSerial_ && comm_->communicator().size() > 1;
460 void prepareFlexibleSolver()
462 OPM_TIMEBLOCK(flexibleSolverPrepare);
465 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
466 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
468 std::function<Vector()> weightCalculator = this->getWeightsCalculator(prm_[activeSolverNum_], getMatrix(), pressureIndex);
469 OPM_TIMEBLOCK(flexibleSolverCreate);
470 flexibleSolver_[activeSolverNum_].create(getMatrix(),
472 prm_[activeSolverNum_],
480 OPM_TIMEBLOCK(flexibleSolverUpdate);
481 flexibleSolver_[activeSolverNum_].pre_->update();
492 if (flexibleSolver_.empty()) {
495 if (!flexibleSolver_[activeSolverNum_].solver_) {
498 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 0) {
502 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 1) {
504 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
505 return newton_iteration == 0;
507 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 2) {
511 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 3) {
515 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 4) {
517 const int step = this->parameters_[activeSolverNum_].cpr_reuse_interval_;
518 const bool create = ((solveCount_ % step) == 0);
523 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
524 std::string msg =
"Invalid value: " + std::to_string(this->parameters_[activeSolverNum_].cpr_reuse_setup_)
525 +
" for --cpr-reuse-setup parameter, run with --help to see allowed values.";
529 throw std::runtime_error(msg);
539 std::function<Vector()> getWeightsCalculator(
const PropertyTree& prm,
540 const Matrix& matrix,
541 std::size_t pressIndex)
const
543 std::function<Vector()> weightsCalculator;
545 using namespace std::string_literals;
547 auto preconditionerType = prm.get(
"preconditioner.type"s,
"cpr"s);
548 if (preconditionerType ==
"cpr" || preconditionerType ==
"cprt"
549 || preconditionerType ==
"cprw" || preconditionerType ==
"cprwt") {
550 const bool transpose = preconditionerType ==
"cprt" || preconditionerType ==
"cprwt";
551 const auto weightsType = prm.get(
"preconditioner.weight_type"s,
"quasiimpes"s);
552 if (weightsType ==
"quasiimpes") {
556 weightsCalculator = [matrix, transpose, pressIndex]() {
557 return Amg::getQuasiImpesWeights<Matrix, Vector>(matrix,
561 }
else if ( weightsType ==
"trueimpes" ) {
565 Vector weights(rhs_->size());
566 ElementContext elemCtx(simulator_);
567 Amg::getTrueImpesWeights(pressIndex, weights,
568 simulator_.vanguard().gridView(),
569 elemCtx, simulator_.model(),
570 ThreadManager::threadId());
573 }
else if (weightsType ==
"trueimpesanalytic" ) {
577 Vector weights(rhs_->size());
578 ElementContext elemCtx(simulator_);
579 Amg::getTrueImpesWeightsAnalytic(pressIndex, weights,
580 simulator_.vanguard().gridView(),
581 elemCtx, simulator_.model(),
582 ThreadManager::threadId());
586 OPM_THROW(std::invalid_argument,
587 "Weights type " + weightsType +
588 "not implemented for cpr."
589 " Please use quasiimpes, trueimpes or trueimpesanalytic.");
592 return weightsCalculator;
601 const Matrix& getMatrix()
const
606 const Simulator& simulator_;
607 mutable int iterations_;
608 mutable int solveCount_;
609 mutable bool converged_;
610 std::any parallelInformation_;
616 int activeSolverNum_ = 0;
617 std::vector<detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType>> flexibleSolver_;
618 std::vector<int> overlapRows_;
619 std::vector<int> interiorRows_;
623 std::vector<FlowLinearSolverParameters> parameters_;
624 bool forceSerial_ =
false;
625 std::vector<PropertyTree> prm_;
627 std::shared_ptr< CommunicationType > comm_;