209 using Matrix =
typename SparseMatrixAdapter::IstlMatrix;
212 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
213 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
220 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
222 using CommunicationType = Dune::CollectiveCommunication<int>;
226 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
228 static void registerParameters()
230 FlowLinearSolverParameters::registerParameters<TypeTag>();
246 comm_.reset(
new CommunicationType( simulator_.vanguard().grid().comm() ) );
248 parameters_.template
init<TypeTag>(simulator_.vanguard().eclState().getSimulationConfig().useCPR());
253#if COMPILE_BDA_BRIDGE
255 std::string accelerator_mode =
EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
256 if ((simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode !=
"none")) {
258 OpmLog::warning(
"Cannot use GPU with MPI, GPU are disabled");
260 accelerator_mode =
"none";
262 const int platformID =
EWOMS_GET_PARAM(TypeTag,
int, OpenclPlatformId);
265 const double tolerance =
EWOMS_GET_PARAM(TypeTag,
double, LinearSolverReduction);
266 const bool opencl_ilu_parallel =
EWOMS_GET_PARAM(TypeTag,
bool, OpenclIluParallel);
269 bdaBridge = std::make_unique<detail::BdaSolverInfo<Matrix,Vector>>(accelerator_mode,
279 if (
EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) !=
"none") {
280 OPM_THROW(std::logic_error,
"Cannot use accelerated solver since CUDA, OpenCL and amgcl were not found by cmake");
289 detail::findOverlapAndInterior(simulator_.vanguard().grid(),
elemMapper, overlapRows_, interiorRows_);
290 useWellConn_ =
EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
293 const std::string
msg =
"The linear solver no longer supports --owner-cells-first=false.";
300 flexibleSolver_.interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(),
true);
304 std::ostringstream
os;
305 os <<
"Property tree for linear solver:\n";
306 prm_.write_json(
os,
true);
307 OpmLog::note(
os.str());
316 void prepare(
const SparseMatrixAdapter& M, Vector& b)
318 OPM_TIMEBLOCK(istlSolverEbosPrepare);
319 static bool firstcall =
true;
322 const size_t size = M.istlMatrix().N();
323 detail::copyParValues(parallelInformation_, size, *comm_);
332 matrix_ =
const_cast<Matrix*
>(&M.istlMatrix());
334 useWellConn_ = EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
337 bdaBridge->numJacobiBlocks_ = EWOMS_GET_PARAM(TypeTag,
int, NumJacobiBlocks);
338 bdaBridge->prepare(simulator_.vanguard().grid(),
339 simulator_.vanguard().cartesianIndexMapper(),
340 simulator_.vanguard().schedule().getWellsatEnd(),
341 simulator_.vanguard().cellPartition(),
342 getMatrix().nonzeroes(), useWellConn_);
346 if ( &(M.istlMatrix()) != matrix_ ) {
347 OPM_THROW(std::logic_error,
348 "Matrix objects are expected to be reused when reassembling!");
353 if (isParallel() && prm_.get<std::string>(
"preconditioner.type") !=
"ParOverILU0") {
354 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
356#if COMPILE_BDA_BRIDGE
357 if(!bdaBridge->gpuActive()){
358 prepareFlexibleSolver();
361 prepareFlexibleSolver();
368 void setResidual(Vector& )
373 void getResidual(Vector& b)
const
378 void setMatrix(
const SparseMatrixAdapter& )
383 bool solve(Vector& x)
385 OPM_TIMEBLOCK(istlSolverEbosSolve);
388 const int verbosity = prm_.get<
int>(
"verbosity", 0);
389 const bool write_matrix = verbosity > 10;
391 Helper::writeSystem(simulator_,
398 Dune::InverseOperatorResult result;
400#if COMPILE_BDA_BRIDGE
401 std::function<void(WellContributions&)> getContribs =
402 [
this](WellContributions& w)
404 this->simulator_.problem().wellModel().getWellContributions(w);
406 if (!bdaBridge->apply(*rhs_, useWellConn_, getContribs,
407 simulator_.gridView().comm().rank(),
408 const_cast<Matrix&
>(this->getMatrix()),
412 OPM_TIMEBLOCK(flexibleSolverApply);
413#if COMPILE_BDA_BRIDGE
414 if(bdaBridge->gpuActive()){
415 prepareFlexibleSolver();
418 assert(flexibleSolver_.solver_);
419 flexibleSolver_.solver_->apply(x, *rhs_, result);
423 checkConvergence(result);
443 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
446 void checkConvergence(
const Dune::InverseOperatorResult&
result )
const
449 iterations_ =
result.iterations;
450 converged_ =
result.converged;
452 if(
result.reduction < parameters_.relaxed_linear_solver_reduction_){
453 std::stringstream
ss;
454 ss<<
"Full linear solver tolerance not achieved. The reduction is:" <<
result.reduction
455 <<
" after " <<
result.iterations <<
" iterations ";
456 OpmLog::warning(
ss.str());
461 if (!parameters_.ignoreConvergenceFailure_ && !converged_) {
462 const std::string msg(
"Convergence failure for linear solver.");
463 OPM_THROW_NOLOG(NumericalProblem, msg);
468 bool isParallel()
const {
470 return comm_->communicator().size() > 1;
476 void prepareFlexibleSolver()
478 OPM_TIMEBLOCK(flexibleSolverPrepare);
480 std::function<Vector()> trueFunc =
483 return this->getTrueImpesWeights(pressureIndex);
487 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
488 flexibleSolver_.wellOperator_ = std::move(wellOp);
490 OPM_TIMEBLOCK(flexibleSolverCreate);
491 flexibleSolver_.create(getMatrix(),
500 OPM_TIMEBLOCK(flexibleSolverUpdate);
501 flexibleSolver_.pre_->update();
512 if (!flexibleSolver_.solver_) {
515 if (this->parameters_.cpr_reuse_setup_ == 0) {
519 if (this->parameters_.cpr_reuse_setup_ == 1) {
521 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
524 if (this->parameters_.cpr_reuse_setup_ == 2) {
528 if (this->parameters_.cpr_reuse_setup_ == 3) {
532 if (this->parameters_.cpr_reuse_setup_ == 4) {
534 const int step = this->parameters_.cpr_reuse_interval_;
535 const bool create = ((calls_ % step) == 0);
540 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
541 std::string
msg =
"Invalid value: " + std::to_string(this->parameters_.cpr_reuse_setup_)
542 +
" for --cpr-reuse-setup parameter, run with --help to see allowed values.";
546 throw std::runtime_error(
msg);
556 Vector getTrueImpesWeights(
int pressureVarIndex)
const
559 Vector weights(rhs_->size());
560 ElementContext
elemCtx(simulator_);
561 Amg::getTrueImpesWeights(pressureVarIndex, weights,
562 simulator_.vanguard().gridView(),
564 ThreadManager::threadId());
573 const Matrix& getMatrix()
const
578 const Simulator& simulator_;
579 mutable int iterations_;
581 mutable bool converged_;
582 std::any parallelInformation_;
588#if COMPILE_BDA_BRIDGE
589 std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge;
592 detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType> flexibleSolver_;
593 std::vector<int> overlapRows_;
594 std::vector<int> interiorRows_;
598 FlowLinearSolverParameters parameters_;
601 std::shared_ptr< CommunicationType > comm_;