My Project
Loading...
Searching...
No Matches
ISTLSolverBda.hpp
1/*
2 Copyright 2016 IRIS AS
3 Copyright 2019, 2020 Equinor ASA
4 Copyright 2020 SINTEF Digital, Mathematics and Cybernetics
5
6 This file is part of the Open Porous Media project (OPM).
7
8 OPM is free software: you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12
13 OPM is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with OPM. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#ifndef OPM_ISTLSOLVER_WITH_BDA_HEADER_INCLUDED
23#define OPM_ISTLSOLVER_WITH_BDA_HEADER_INCLUDED
24
25#include <opm/simulators/linalg/ISTLSolver.hpp>
26
27#include <cstddef>
28#include <memory>
29#include <set>
30#include <string>
31#include <vector>
32
33namespace Opm {
34
35class Well;
36
37template<class Matrix, class Vector, int block_size> class BdaBridge;
38class WellContributions;
39namespace detail {
40
41template<class Matrix, class Vector>
43{
44 using WellContribFunc = std::function<void(WellContributions&)>;
46
47 BdaSolverInfo(const std::string& accelerator_mode,
48 const int linear_solver_verbosity,
49 const int maxit,
50 const double tolerance,
51 const int platformID,
52 const int deviceID,
53 const bool opencl_ilu_parallel,
54 const std::string& linsolver);
55
57
58 template<class Grid>
59 void prepare(const Grid& grid,
60 const Dune::CartesianIndexMapper<Grid>& cartMapper,
61 const std::vector<Well>& wellsForConn,
62 const std::vector<int>& cellPartition,
63 const std::size_t nonzeroes,
64 const bool useWellConn);
65
66 bool apply(Vector& rhs,
67 const bool useWellConn,
68 WellContribFunc getContribs,
69 const int rank,
70 Matrix& matrix,
71 Vector& x,
72 Dune::InverseOperatorResult& result);
73
74 bool gpuActive();
75
76 int numJacobiBlocks_ = 0;
77
78private:
81 template<class Grid>
82 void blockJacobiAdjacency(const Grid& grid,
83 const std::vector<int>& cell_part,
84 std::size_t nonzeroes);
85
86 void copyMatToBlockJac(const Matrix& mat, Matrix& blockJac);
87
88 std::unique_ptr<Bridge> bridge_;
89 std::string accelerator_mode_;
90 std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
91 std::vector<std::set<int>> wellConnectionsGraph_;
92};
93
94}
95
100template <class TypeTag>
101class ISTLSolverBda : public ISTLSolver<TypeTag>
102{
103protected:
104 using ParentType = ISTLSolver<TypeTag>;
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;
121
122#if HAVE_MPI
123 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
124#else
125 using CommunicationType = Dune::CollectiveCommunication<int>;
126#endif
127
128public:
129 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
130
135 ISTLSolverBda(const Simulator& simulator, const FlowLinearSolverParameters& parameters)
136 : ParentType(simulator, parameters)
137 {
138 initializeBda();
139 }
140
143 explicit ISTLSolverBda(const Simulator& simulator)
144 : ParentType(simulator)
145 {
146 initializeBda();
147 }
148
149 void initializeBda()
150 {
151 OPM_TIMEBLOCK(initializeBda);
152
153 std::string accelerator_mode = Parameters::get<TypeTag, Properties::AcceleratorMode>();
154 // Force accelerator mode to none if using MPI.
155 if ((this->simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode != "none")) {
156 const bool on_io_rank = (this->simulator_.gridView().comm().rank() == 0);
157 if (on_io_rank) {
158 OpmLog::warning("Cannot use AcceleratorMode feature with MPI, setting AcceleratorMode to 'none'.");
159 }
160 accelerator_mode = "none";
161 }
162
163 if (accelerator_mode == "none") {
164 return;
165 }
166
167 // Initialize the BdaBridge
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,
177 maxit,
178 tolerance,
179 platformID,
180 deviceID,
181 opencl_ilu_parallel,
182 linsolver);
183 }
184
185 void prepare(const Matrix& M, Vector& b)
186 {
187 OPM_TIMEBLOCK(prepare);
188 [[maybe_unused]] const bool firstcall = (this->matrix_ == nullptr);
189
190 // Avoid performing the decomposition on CPU when we also do it on GPU,
191 // but we do need to initialize the pointers.
192 if (bdaBridge_) {
193 ParentType::initPrepare(M,b);
194 } else {
195 ParentType::prepare(M,b);
196 }
197
198#if HAVE_OPENCL || HAVE_ROCSPARSE || HAVE_CUDA
199 // update matrix entries for solvers.
200 if (firstcall && bdaBridge_) {
201 // model will not change the matrix object. Hence simply store a pointer
202 // to the original one with a deleter that does nothing.
203 // Outch! We need to be able to scale the linear system! Hence const_cast
204 // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver)
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_);
211 }
212#endif
213 }
214
215
216 void setResidual(Vector& /* b */)
217 {
218 // rhs_ = &b; // Must be handled in prepare() instead.
219 }
220
221 void getResidual(Vector& b) const
222 {
223 b = *(this->rhs_);
224 }
225
226 void setMatrix(const SparseMatrixAdapter& /* M */)
227 {
228 // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
229 }
230
231 bool solve(Vector& x)
232 {
233 if (!bdaBridge_) {
234 return ParentType::solve(x);
235 }
236
237 OPM_TIMEBLOCK(istlSolverBdaSolve);
238 this->solveCount_ += 1;
239 // Write linear system if asked for.
240 const int verbosity = this->prm_[this->activeSolverNum_].template get<int>("verbosity", 0);
241 const bool write_matrix = verbosity > 10;
242 if (write_matrix) {
243 Helper::writeSystem(this->simulator_, //simulator is only used to get names
244 this->getMatrix(),
245 *(this->rhs_),
246 this->comm_.get());
247 }
248
249 // Solve system.
250 Dune::InverseOperatorResult result;
251
252 std::function<void(WellContributions&)> getContribs =
253 [this](WellContributions& w)
254 {
255 this->simulator_.problem().wellModel().getWellContributions(w);
256 };
257 if (!bdaBridge_->apply(*(this->rhs_), this->useWellConn_, getContribs,
258 this->simulator_.gridView().comm().rank(),
259 const_cast<Matrix&>(this->getMatrix()),
260 x, result))
261 {
262 if(bdaBridge_->gpuActive()){
263 // bda solve fails use istl solver setup need to be done since it is not setup in prepare
264 ParentType::prepareFlexibleSolver();
265 }
266 assert(this->flexibleSolver_[this->activeSolverNum_].solver_);
267 this->flexibleSolver_[this->activeSolverNum_].solver_->apply(x, *(this->rhs_), result);
268 }
269
270 // Check convergence, iterations etc.
271 this->checkConvergence(result);
272
273 return this->converged_;
274 }
275
276protected:
277 std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge_;
278}; // end ISTLSolver
279
280} // namespace Opm
281
282#endif // OPM_ISTLSOLVER_WITH_BDA_HEADER_INCLUDED
Definition CollectDataOnIORank.hpp:49
Interface class adding the update() method to the preconditioner interface.
Definition PreconditionerWithUpdate.hpp:32
BdaBridge acts as interface between opm-simulators with the BdaSolvers.
Definition BdaBridge.hpp:37
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition ISTLSolverBda.hpp:102
ISTLSolverBda(const Simulator &simulator, const FlowLinearSolverParameters &parameters)
Construct a system solver.
Definition ISTLSolverBda.hpp:135
ISTLSolverBda(const Simulator &simulator)
Construct a system solver.
Definition ISTLSolverBda.hpp:143
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition ISTLSolver.hpp:143
This class serves to eliminate the need to include the WellContributions into the matrix (with –matri...
Definition WellContributions.hpp:52
Definition WellOperators.hpp:67
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
This class carries all parameters for the NewtonIterationBlackoilInterleaved class.
Definition FlowLinearSolverParameters.hpp:237
Definition ISTLSolverBda.hpp:43