48 const int linear_solver_verbosity,
50 const double tolerance,
53 const bool opencl_ilu_parallel,
54 const std::string& linsolver);
59 void prepare(
const Grid& grid,
61 const std::vector<Well>& wellsForConn,
62 const std::vector<int>& cellPartition,
63 const std::size_t nonzeroes,
64 const bool useWellConn);
66 bool apply(Vector& rhs,
67 const bool useWellConn,
68 WellContribFunc getContribs,
72 Dune::InverseOperatorResult& result);
76 int numJacobiBlocks_ = 0;
82 void blockJacobiAdjacency(
const Grid& grid,
83 const std::vector<int>& cell_part,
84 std::size_t nonzeroes);
86 void copyMatToBlockJac(
const Matrix& mat, Matrix& blockJac);
88 std::unique_ptr<Bridge> bridge_;
89 std::string accelerator_mode_;
90 std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
91 std::vector<std::set<int>> wellConnectionsGraph_;
105 using GridView = GetPropType<TypeTag, Properties::GridView>;
106 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
107 using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
108 using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
109 using Indices = GetPropType<TypeTag, Properties::Indices>;
110 using WellModel = GetPropType<TypeTag, Properties::WellModel>;
111 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
112 using Matrix =
typename SparseMatrixAdapter::IstlMatrix;
113 using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
114 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
115 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
116 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
119 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
120 constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
123 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
125 using CommunicationType = Dune::CollectiveCommunication<int>;
129 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
136 : ParentType(simulator, parameters)
144 : ParentType(simulator)
151 OPM_TIMEBLOCK(initializeBda);
153 std::string accelerator_mode = Parameters::get<TypeTag, Properties::AcceleratorMode>();
155 if ((this->simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode !=
"none")) {
156 const bool on_io_rank = (this->simulator_.gridView().comm().rank() == 0);
158 OpmLog::warning(
"Cannot use AcceleratorMode feature with MPI, setting AcceleratorMode to 'none'.");
160 accelerator_mode =
"none";
163 if (accelerator_mode ==
"none") {
168 const int platformID = Parameters::get<TypeTag, Properties::OpenclPlatformId>();
169 const int deviceID = Parameters::get<TypeTag, Properties::BdaDeviceId>();
170 const int maxit = Parameters::get<TypeTag, Properties::LinearSolverMaxIter>();
171 const double tolerance = Parameters::get<TypeTag, Properties::LinearSolverReduction>();
172 const bool opencl_ilu_parallel = Parameters::get<TypeTag, Properties::OpenclIluParallel>();
173 const int linear_solver_verbosity = this->parameters_[0].linear_solver_verbosity_;
174 std::string linsolver = Parameters::get<TypeTag, Properties::LinearSolver>();
175 bdaBridge_ = std::make_unique<detail::BdaSolverInfo<Matrix,Vector>>(accelerator_mode,
176 linear_solver_verbosity,
185 void prepare(
const Matrix& M, Vector& b)
187 OPM_TIMEBLOCK(prepare);
188 [[maybe_unused]]
const bool firstcall = (this->matrix_ ==
nullptr);
193 ParentType::initPrepare(M,b);
195 ParentType::prepare(M,b);
198#if HAVE_OPENCL || HAVE_ROCSPARSE || HAVE_CUDA
200 if (firstcall && bdaBridge_) {
205 bdaBridge_->numJacobiBlocks_ = Parameters::get<TypeTag, Properties::NumJacobiBlocks>();
206 bdaBridge_->prepare(this->simulator_.vanguard().grid(),
207 this->simulator_.vanguard().cartesianIndexMapper(),
208 this->simulator_.vanguard().schedule().getWellsatEnd(),
209 this->simulator_.vanguard().cellPartition(),
210 this->getMatrix().nonzeroes(), this->useWellConn_);
216 void setResidual(Vector& )
221 void getResidual(Vector& b)
const
226 void setMatrix(
const SparseMatrixAdapter& )
231 bool solve(Vector& x)
234 return ParentType::solve(x);
237 OPM_TIMEBLOCK(istlSolverBdaSolve);
238 this->solveCount_ += 1;
240 const int verbosity = this->prm_[this->activeSolverNum_].template get<int>(
"verbosity", 0);
241 const bool write_matrix = verbosity > 10;
243 Helper::writeSystem(this->simulator_,
250 Dune::InverseOperatorResult result;
252 std::function<void(WellContributions&)> getContribs =
253 [
this](WellContributions& w)
255 this->simulator_.problem().wellModel().getWellContributions(w);
257 if (!bdaBridge_->apply(*(this->rhs_), this->useWellConn_, getContribs,
258 this->simulator_.gridView().comm().rank(),
259 const_cast<Matrix&
>(this->getMatrix()),
262 if(bdaBridge_->gpuActive()){
264 ParentType::prepareFlexibleSolver();
266 assert(this->flexibleSolver_[this->activeSolverNum_].solver_);
267 this->flexibleSolver_[this->activeSolverNum_].solver_->apply(x, *(this->rhs_), result);
271 this->checkConvergence(result);
273 return this->converged_;
277 std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge_;