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);
281 std::cerr<<
"debug: ["<<i<<
"] : " << os.str() <<std::endl;
283 OpmLog::note(os.str());
292 void setActiveSolver(
const int num)
294 if (num >
static_cast<int>(prm_.size()) - 1) {
295 OPM_THROW(std::logic_error,
"Solver number " + std::to_string(num) +
" not available.");
297 activeSolverNum_ = num;
298 auto cc = Dune::MPIHelper::getCollectiveCommunication();
299 if (cc.rank() == 0) {
300 OpmLog::debug(
"Active solver = " + std::to_string(activeSolverNum_)
301 +
" (" + parameters_[activeSolverNum_].linsolver_ +
")");
305 int numAvailableSolvers()
307 return flexibleSolver_.size();
310 void prepare(
const SparseMatrixAdapter& M, Vector& b)
312 prepare(M.istlMatrix(), b);
315 void prepare(
const Matrix& M, Vector& b)
317 OPM_TIMEBLOCK(istlSolverEbosPrepare);
318 const bool firstcall = (matrix_ ==
nullptr);
320 if (firstcall && isParallel()) {
321 const std::size_t size = M.N();
322 detail::copyParValues(parallelInformation_, size, *comm_);
331 matrix_ =
const_cast<Matrix*
>(&M);
333 useWellConn_ = EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
337 if ( &M != matrix_ ) {
338 OPM_THROW(std::logic_error,
339 "Matrix objects are expected to be reused when reassembling!");
345 if (isParallel() && prm_[activeSolverNum_].
template get<std::string>(
"preconditioner.type") !=
"ParOverILU0") {
346 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
348 prepareFlexibleSolver();
352 void setResidual(Vector& )
357 void getResidual(Vector& b)
const
362 void setMatrix(
const SparseMatrixAdapter& )
367 int getSolveCount()
const {
371 void resetSolveCount() {
375 bool solve(Vector& x)
377 OPM_TIMEBLOCK(istlSolverEbosSolve);
380 const int verbosity = prm_[activeSolverNum_].get(
"verbosity", 0);
381 const bool write_matrix = verbosity > 10;
383 Helper::writeSystem(simulator_,
390 Dune::InverseOperatorResult result;
392 OPM_TIMEBLOCK(flexibleSolverApply);
393 assert(flexibleSolver_[activeSolverNum_].solver_);
394 flexibleSolver_[activeSolverNum_].solver_->apply(x, *rhs_, result);
398 checkConvergence(result);
416 const CommunicationType* comm()
const {
return comm_.get(); }
420 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
423 void checkConvergence(
const Dune::InverseOperatorResult& result )
const
426 iterations_ = result.iterations;
427 converged_ = result.converged;
429 if(result.reduction < parameters_[activeSolverNum_].relaxed_linear_solver_reduction_){
430 std::stringstream ss;
431 ss<<
"Full linear solver tolerance not achieved. The reduction is:" << result.reduction
432 <<
" after " << result.iterations <<
" iterations ";
433 OpmLog::warning(ss.str());
438 if (!parameters_[activeSolverNum_].ignoreConvergenceFailure_ && !converged_) {
439 const std::string msg(
"Convergence failure for linear solver.");
440 OPM_THROW_NOLOG(NumericalProblem, msg);
445 bool isParallel()
const {
447 return !forceSerial_ && comm_->communicator().size() > 1;
453 void prepareFlexibleSolver()
455 OPM_TIMEBLOCK(flexibleSolverPrepare);
457 std::function<Vector()> trueFunc =
460 return this->getTrueImpesWeights(pressureIndex);
464 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
465 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
467 OPM_TIMEBLOCK(flexibleSolverCreate);
468 flexibleSolver_[activeSolverNum_].create(getMatrix(),
470 prm_[activeSolverNum_],
478 OPM_TIMEBLOCK(flexibleSolverUpdate);
479 flexibleSolver_[activeSolverNum_].pre_->update();
490 if (flexibleSolver_.empty()) {
493 if (!flexibleSolver_[activeSolverNum_].solver_) {
496 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 0) {
500 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 1) {
502 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
505 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 2) {
509 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 3) {
513 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 4) {
515 const int step = this->parameters_[activeSolverNum_].cpr_reuse_interval_;
516 const bool create = ((solveCount_ % step) == 0);
521 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
522 std::string
msg =
"Invalid value: " + std::to_string(this->parameters_[activeSolverNum_].cpr_reuse_setup_)
523 +
" for --cpr-reuse-setup parameter, run with --help to see allowed values.";
527 throw std::runtime_error(
msg);
537 Vector getTrueImpesWeights(
int pressureVarIndex)
const
540 Vector weights(rhs_->size());
541 ElementContext
elemCtx(simulator_);
542 Amg::getTrueImpesWeights(pressureVarIndex, weights,
543 simulator_.vanguard().gridView(),
545 ThreadManager::threadId());
554 const Matrix& getMatrix()
const
559 const Simulator& simulator_;
560 mutable int iterations_;
561 mutable int solveCount_;
562 mutable bool converged_;
563 std::any parallelInformation_;
569 int activeSolverNum_ = 0;
570 std::vector<detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType>> flexibleSolver_;
571 std::vector<int> overlapRows_;
572 std::vector<int> interiorRows_;
576 std::vector<FlowLinearSolverParameters> parameters_;
577 bool forceSerial_ =
false;
578 std::vector<PropertyTree> prm_;
580 std::shared_ptr< CommunicationType > comm_;