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 auto cc = Dune::MPIHelper::getCollectiveCommunication();
298 if (cc.rank() == 0) {
299 OpmLog::debug(
"Active solver = " + std::to_string(activeSolverNum_)
300 +
" (" + parameters_[activeSolverNum_].linsolver_ +
")");
304 int numAvailableSolvers()
306 return flexibleSolver_.size();
309 void prepare(
const SparseMatrixAdapter& M, Vector& b)
311 prepare(M.istlMatrix(), b);
314 void prepare(
const Matrix& M, Vector& b)
316 OPM_TIMEBLOCK(istlSolverEbosPrepare);
317 const bool firstcall = (matrix_ ==
nullptr);
319 if (firstcall && isParallel()) {
320 const std::size_t size = M.N();
321 detail::copyParValues(parallelInformation_, size, *comm_);
330 matrix_ =
const_cast<Matrix*
>(&M);
332 useWellConn_ = EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
336 if ( &M != matrix_ ) {
337 OPM_THROW(std::logic_error,
338 "Matrix objects are expected to be reused when reassembling!");
344 if (isParallel() && prm_[activeSolverNum_].
template get<std::string>(
"preconditioner.type") !=
"ParOverILU0") {
345 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
347 prepareFlexibleSolver();
351 void setResidual(Vector& )
356 void getResidual(Vector& b)
const
361 void setMatrix(
const SparseMatrixAdapter& )
366 int getSolveCount()
const {
370 void resetSolveCount() {
374 bool solve(Vector& x)
376 OPM_TIMEBLOCK(istlSolverEbosSolve);
379 const int verbosity = prm_[activeSolverNum_].get(
"verbosity", 0);
380 const bool write_matrix = verbosity > 10;
382 Helper::writeSystem(simulator_,
389 Dune::InverseOperatorResult result;
391 OPM_TIMEBLOCK(flexibleSolverApply);
392 assert(flexibleSolver_[activeSolverNum_].solver_);
393 flexibleSolver_[activeSolverNum_].solver_->apply(x, *rhs_, result);
397 checkConvergence(result);
415 const CommunicationType* comm()
const {
return comm_.get(); }
419 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
422 void checkConvergence(
const Dune::InverseOperatorResult& result )
const
425 iterations_ = result.iterations;
426 converged_ = result.converged;
428 if(result.reduction < parameters_[activeSolverNum_].relaxed_linear_solver_reduction_){
429 std::stringstream ss;
430 ss<<
"Full linear solver tolerance not achieved. The reduction is:" << result.reduction
431 <<
" after " << result.iterations <<
" iterations ";
432 OpmLog::warning(ss.str());
437 if (!parameters_[activeSolverNum_].ignoreConvergenceFailure_ && !converged_) {
438 const std::string msg(
"Convergence failure for linear solver.");
439 OPM_THROW_NOLOG(NumericalProblem, msg);
444 bool isParallel()
const {
446 return !forceSerial_ && comm_->communicator().size() > 1;
452 void prepareFlexibleSolver()
454 OPM_TIMEBLOCK(flexibleSolverPrepare);
456 std::function<Vector()> trueFunc =
459 return this->getTrueImpesWeights(pressureIndex);
463 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
464 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
466 OPM_TIMEBLOCK(flexibleSolverCreate);
467 flexibleSolver_[activeSolverNum_].create(getMatrix(),
469 prm_[activeSolverNum_],
477 OPM_TIMEBLOCK(flexibleSolverUpdate);
478 flexibleSolver_[activeSolverNum_].pre_->update();
489 if (flexibleSolver_.empty()) {
492 if (!flexibleSolver_[activeSolverNum_].solver_) {
495 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 0) {
499 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 1) {
501 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
504 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 2) {
508 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 3) {
512 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 4) {
514 const int step = this->parameters_[activeSolverNum_].cpr_reuse_interval_;
515 const bool create = ((solveCount_ % step) == 0);
520 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
521 std::string
msg =
"Invalid value: " + std::to_string(this->parameters_[activeSolverNum_].cpr_reuse_setup_)
522 +
" for --cpr-reuse-setup parameter, run with --help to see allowed values.";
526 throw std::runtime_error(
msg);
536 Vector getTrueImpesWeights(
int pressureVarIndex)
const
539 Vector weights(rhs_->size());
540 ElementContext
elemCtx(simulator_);
541 Amg::getTrueImpesWeights(pressureVarIndex, weights,
542 simulator_.vanguard().gridView(),
544 ThreadManager::threadId());
553 const Matrix& getMatrix()
const
558 const Simulator& simulator_;
559 mutable int iterations_;
560 mutable int solveCount_;
561 mutable bool converged_;
562 std::any parallelInformation_;
568 int activeSolverNum_ = 0;
569 std::vector<detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType>> flexibleSolver_;
570 std::vector<int> overlapRows_;
571 std::vector<int> interiorRows_;
575 std::vector<FlowLinearSolverParameters> parameters_;
576 bool forceSerial_ =
false;
577 std::vector<PropertyTree> prm_;
579 std::shared_ptr< CommunicationType > comm_;