153 using Matrix =
typename SparseMatrixAdapter::IstlMatrix;
156 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
157 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
164 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
166 using CommunicationType = Dune::CollectiveCommunication<int>;
170 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
172 static void registerParameters()
174 FlowLinearSolverParameters::registerParameters<TypeTag>();
204 parameters_.resize(1);
205 parameters_[0].template
init<TypeTag>(simulator_.vanguard().eclState().getSimulationConfig().useCPR());
213 if (parameters_[0].linsolver_ ==
"hybrid") {
222 para.init<TypeTag>(
false);
223 para.linsolver_ =
"cprw";
224 parameters_.push_back(
para);
230 FlowLinearSolverParameters para;
231 para.init<TypeTag>(
false);
232 para.linsolver_ =
"ilu0";
233 parameters_.push_back(para);
235 EWOMS_PARAM_IS_SET(TypeTag,
int, LinearSolverMaxIter),
236 EWOMS_PARAM_IS_SET(TypeTag,
double, LinearSolverReduction)));
241 assert(parameters_.size() == 1);
242 assert(prm_.empty());
244 EWOMS_PARAM_IS_SET(TypeTag,
int, LinearSolverMaxIter),
245 EWOMS_PARAM_IS_SET(TypeTag,
double, LinearSolverReduction)));
247 flexibleSolver_.resize(prm_.size());
249 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
251 comm_.reset(
new CommunicationType( simulator_.vanguard().grid().comm() ) );
253 extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
258 ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
259 detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
260 useWellConn_ = EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
261 const bool ownersFirst = EWOMS_GET_PARAM(TypeTag,
bool, OwnerCellsFirst);
263 const std::string msg =
"The linear solver no longer supports --owner-cells-first=false.";
267 OPM_THROW_NOLOG(std::runtime_error, msg);
270 const int interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(),
true);
271 for (
auto& f : flexibleSolver_) {
272 f.interiorCellNum_ = interiorCellNum_;
276 if (on_io_rank && parameters_[activeSolverNum_].linear_solver_print_json_definition_) {
277 std::ostringstream os;
278 os <<
"Property tree for linear solvers:\n";
279 for (std::size_t i = 0; i<prm_.size(); i++) {
280 prm_[i].write_json(os,
true);
282 OpmLog::note(os.str());
291 void setActiveSolver(
const int num)
293 if (num >
static_cast<int>(prm_.size()) - 1) {
294 OPM_THROW(std::logic_error,
"Solver number " + std::to_string(num) +
" not available.");
296 activeSolverNum_ = num;
297 if (simulator_.gridView().comm().rank() == 0) {
298 OpmLog::debug(
"Active solver = " + std::to_string(activeSolverNum_)
299 +
" (" + parameters_[activeSolverNum_].linsolver_ +
")");
303 int numAvailableSolvers()
305 return flexibleSolver_.size();
308 void prepare(
const SparseMatrixAdapter& M, Vector& b)
310 prepare(M.istlMatrix(), b);
313 void prepare(
const Matrix& M, Vector& b)
315 OPM_TIMEBLOCK(istlSolverEbosPrepare);
316 const bool firstcall = (matrix_ ==
nullptr);
318 if (firstcall && isParallel()) {
319 const std::size_t size = M.N();
320 detail::copyParValues(parallelInformation_, size, *comm_);
329 matrix_ =
const_cast<Matrix*
>(&M);
331 useWellConn_ = EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
335 if ( &M != matrix_ ) {
336 OPM_THROW(std::logic_error,
337 "Matrix objects are expected to be reused when reassembling!");
343 if (isParallel() && prm_[activeSolverNum_].
template get<std::string>(
"preconditioner.type") !=
"ParOverILU0") {
344 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
346 prepareFlexibleSolver();
350 void setResidual(Vector& )
355 void getResidual(Vector& b)
const
360 void setMatrix(
const SparseMatrixAdapter& )
365 int getSolveCount()
const {
369 void resetSolveCount() {
373 bool solve(Vector& x)
375 OPM_TIMEBLOCK(istlSolverEbosSolve);
378 const int verbosity = prm_[activeSolverNum_].get(
"verbosity", 0);
379 const bool write_matrix = verbosity > 10;
381 Helper::writeSystem(simulator_,
388 Dune::InverseOperatorResult result;
390 OPM_TIMEBLOCK(flexibleSolverApply);
391 assert(flexibleSolver_[activeSolverNum_].solver_);
392 flexibleSolver_[activeSolverNum_].solver_->apply(x, *rhs_, result);
396 checkConvergence(result);
414 const CommunicationType* comm()
const {
return comm_.get(); }
418 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
421 void checkConvergence(
const Dune::InverseOperatorResult& result )
const
424 iterations_ = result.iterations;
425 converged_ = result.converged;
427 if(result.reduction < parameters_[activeSolverNum_].relaxed_linear_solver_reduction_){
428 std::stringstream ss;
429 ss<<
"Full linear solver tolerance not achieved. The reduction is:" << result.reduction
430 <<
" after " << result.iterations <<
" iterations ";
431 OpmLog::warning(ss.str());
436 if (!parameters_[activeSolverNum_].ignoreConvergenceFailure_ && !converged_) {
437 const std::string msg(
"Convergence failure for linear solver.");
438 OPM_THROW_NOLOG(NumericalProblem, msg);
443 bool isParallel()
const {
445 return !forceSerial_ && comm_->communicator().size() > 1;
451 void prepareFlexibleSolver()
453 OPM_TIMEBLOCK(flexibleSolverPrepare);
455 std::function<Vector()> trueFunc =
458 return this->getTrueImpesWeights(pressureIndex);
462 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
463 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
465 OPM_TIMEBLOCK(flexibleSolverCreate);
466 flexibleSolver_[activeSolverNum_].create(getMatrix(),
468 prm_[activeSolverNum_],
476 OPM_TIMEBLOCK(flexibleSolverUpdate);
477 flexibleSolver_[activeSolverNum_].pre_->update();
488 if (flexibleSolver_.empty()) {
491 if (!flexibleSolver_[activeSolverNum_].solver_) {
494 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 0) {
498 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 1) {
500 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
503 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 2) {
507 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 3) {
511 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 4) {
513 const int step = this->parameters_[activeSolverNum_].cpr_reuse_interval_;
514 const bool create = ((solveCount_ % step) == 0);
519 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
520 std::string
msg =
"Invalid value: " + std::to_string(this->parameters_[activeSolverNum_].cpr_reuse_setup_)
521 +
" for --cpr-reuse-setup parameter, run with --help to see allowed values.";
525 throw std::runtime_error(
msg);
535 Vector getTrueImpesWeights(
int pressureVarIndex)
const
538 Vector weights(rhs_->size());
539 ElementContext
elemCtx(simulator_);
540 Amg::getTrueImpesWeights(pressureVarIndex, weights,
541 simulator_.vanguard().gridView(),
543 ThreadManager::threadId());
552 const Matrix& getMatrix()
const
557 const Simulator& simulator_;
558 mutable int iterations_;
559 mutable int solveCount_;
560 mutable bool converged_;
561 std::any parallelInformation_;
567 int activeSolverNum_ = 0;
568 std::vector<detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType>> flexibleSolver_;
569 std::vector<int> overlapRows_;
570 std::vector<int> interiorRows_;
574 std::vector<FlowLinearSolverParameters> parameters_;
575 bool forceSerial_ =
false;
576 std::vector<PropertyTree> prm_;
578 std::shared_ptr< CommunicationType > comm_;