44 using Scalar =
typename Vector::field_type;
49 const int linear_solver_verbosity,
51 const Scalar tolerance,
54 const bool opencl_ilu_parallel,
55 const std::string& linsolver);
60 void prepare(
const Grid& grid,
62 const std::vector<Well>& wellsForConn,
63 const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
64 const std::vector<int>& cellPartition,
65 const std::size_t nonzeroes,
66 const bool useWellConn);
68 bool apply(Vector& rhs,
69 const bool useWellConn,
70 WellContribFunc getContribs,
74 Dune::InverseOperatorResult& result);
78 int numJacobiBlocks_ = 0;
84 void blockJacobiAdjacency(
const Grid& grid,
85 const std::vector<int>& cell_part,
86 std::size_t nonzeroes);
88 void copyMatToBlockJac(
const Matrix& mat, Matrix& blockJac);
90 std::unique_ptr<Bridge> bridge_;
91 std::string accelerator_mode_;
92 std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
93 std::vector<std::set<int>> wellConnectionsGraph_;
114 using Matrix =
typename SparseMatrixAdapter::IstlMatrix;
117 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
118 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
125 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
127 using CommunicationType = Dune::Communication<int>;
131 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
153 OPM_TIMEBLOCK(initializeBda);
155 std::string accelerator_mode = Parameters::Get<Parameters::AcceleratorMode>();
157 if ((this->simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode !=
"none")) {
158 const bool on_io_rank = (this->simulator_.gridView().comm().rank() == 0);
160 OpmLog::warning(
"Cannot use AcceleratorMode feature with MPI, setting AcceleratorMode to 'none'.");
162 accelerator_mode =
"none";
165 if (accelerator_mode ==
"none") {
170 const int platformID = Parameters::Get<Parameters::OpenclPlatformId>();
171 const int deviceID = Parameters::Get<Parameters::BdaDeviceId>();
172 const int maxit = Parameters::Get<Parameters::LinearSolverMaxIter>();
173 const double tolerance = Parameters::Get<Parameters::LinearSolverReduction>();
174 const bool opencl_ilu_parallel = Parameters::Get<Parameters::OpenclIluParallel>();
175 const int linear_solver_verbosity = this->parameters_[0].linear_solver_verbosity_;
176 std::string linsolver = Parameters::Get<Parameters::LinearSolver>();
177 bdaBridge_ = std::make_unique<detail::BdaSolverInfo<Matrix,Vector>>(accelerator_mode,
178 linear_solver_verbosity,
187 void prepare(
const Matrix& M, Vector& b)
189 OPM_TIMEBLOCK(prepare);
190 [[maybe_unused]]
const bool firstcall = (this->matrix_ ==
nullptr);
195 ParentType::initPrepare(M,b);
197 ParentType::prepare(M,b);
200#if HAVE_OPENCL || HAVE_ROCSPARSE || HAVE_CUDA
202 if (firstcall && bdaBridge_) {
207 bdaBridge_->numJacobiBlocks_ = Parameters::Get<Parameters::NumJacobiBlocks>();
208 bdaBridge_->prepare(this->simulator_.vanguard().grid(),
209 this->simulator_.vanguard().cartesianIndexMapper(),
210 this->simulator_.vanguard().schedule().getWellsatEnd(),
211 this->simulator_.vanguard().schedule().getPossibleFutureConnections(),
212 this->simulator_.vanguard().cellPartition(),
213 this->getMatrix().nonzeroes(), this->useWellConn_);
219 void setResidual(Vector& )
224 void getResidual(Vector& b)
const
229 void setMatrix(
const SparseMatrixAdapter& )
234 bool solve(Vector& x)
237 return ParentType::solve(x);
240 OPM_TIMEBLOCK(istlSolverBdaSolve);
241 this->solveCount_ += 1;
243 const int verbosity = this->prm_[this->activeSolverNum_].template get<int>(
"verbosity", 0);
244 const bool write_matrix = verbosity > 10;
246 Helper::writeSystem(this->simulator_,
253 Dune::InverseOperatorResult result;
255 std::function<void(WellContributions<Scalar>&)> getContribs =
256 [
this](WellContributions<Scalar>& w)
258 this->simulator_.problem().wellModel().getWellContributions(w);
260 if (!bdaBridge_->apply(*(this->rhs_), this->useWellConn_, getContribs,
261 this->simulator_.gridView().comm().rank(),
262 const_cast<Matrix&
>(this->getMatrix()),
265 if(bdaBridge_->gpuActive()){
267 ParentType::prepareFlexibleSolver();
269 assert(this->flexibleSolver_[this->activeSolverNum_].solver_);
270 this->flexibleSolver_[this->activeSolverNum_].solver_->apply(x, *(this->rhs_), result);
274 this->checkConvergence(result);
276 return this->converged_;
280 std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge_;
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235