My Project
Loading...
Searching...
No Matches
ISTLSolverEbosBda.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_EBOS_WITH_BDA_INCLUDED
23#define OPM_ISTLSOLVER_EBOS_WITH_BDA_INCLUDED
24
25#include <opm/simulators/linalg/ISTLSolverEbos.hpp>
26
27#include <cstddef>
28
29namespace Opm {
30
31class Well;
32
33template<class Matrix, class Vector, int block_size> class BdaBridge;
34class WellContributions;
35namespace detail {
36
37template<class Matrix, class Vector>
39{
40 using WellContribFunc = std::function<void(WellContributions&)>;
42
43 BdaSolverInfo(const std::string& accelerator_mode,
45 const int maxit,
46 const double tolerance,
47 const int platformID,
48 const int deviceID,
49 const bool opencl_ilu_parallel,
50 const std::string& linsolver);
51
53
54 template<class Grid>
55 void prepare(const Grid& grid,
57 const std::vector<Well>& wellsForConn,
58 const std::vector<int>& cellPartition,
59 const std::size_t nonzeroes,
60 const bool useWellConn);
61
62 bool apply(Vector& rhs,
63 const bool useWellConn,
64 WellContribFunc getContribs,
65 const int rank,
66 Matrix& matrix,
67 Vector& x,
68 Dune::InverseOperatorResult& result);
69
70 bool gpuActive();
71
72 int numJacobiBlocks_ = 0;
73
74private:
77 template<class Grid>
78 void blockJacobiAdjacency(const Grid& grid,
79 const std::vector<int>& cell_part,
80 std::size_t nonzeroes);
81
82 void copyMatToBlockJac(const Matrix& mat, Matrix& blockJac);
83
84 std::unique_ptr<Bridge> bridge_;
85 std::string accelerator_mode_;
86 std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
87 std::vector<std::set<int>> wellConnectionsGraph_;
88};
89
90}
91
96template <class TypeTag>
97class ISTLSolverEbosBda : public ISTLSolverEbos<TypeTag>
98{
99protected:
108 using Matrix = typename SparseMatrixAdapter::IstlMatrix;
111 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
112 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
116 constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
117
118#if HAVE_MPI
119 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
120#else
121 using CommunicationType = Dune::CollectiveCommunication<int>;
122#endif
123
124public:
125 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
126
136
141 {
142 initializeBda();
143 }
144
145 void initializeBda()
146 {
147 OPM_TIMEBLOCK(initializeBda);
148
149 std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
150 // Force accelerator mode to none if using MPI.
151 if ((this->simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode != "none")) {
152 const bool on_io_rank = (this->simulator_.gridView().comm().rank() == 0);
153 if (on_io_rank) {
154 OpmLog::warning("Cannot use AcceleratorMode feature with MPI, setting AcceleratorMode to 'none'.");
155 }
156 accelerator_mode = "none";
157 }
158
159 if (accelerator_mode == "none") {
160 return;
161 }
162
163 // Initialize the BdaBridge
164 const int platformID = EWOMS_GET_PARAM(TypeTag, int, OpenclPlatformId);
165 const int deviceID = EWOMS_GET_PARAM(TypeTag, int, BdaDeviceId);
166 const int maxit = EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIter);
167 const double tolerance = EWOMS_GET_PARAM(TypeTag, double, LinearSolverReduction);
168 const bool opencl_ilu_parallel = EWOMS_GET_PARAM(TypeTag, bool, OpenclIluParallel);
169 const int linear_solver_verbosity = this->parameters_[0].linear_solver_verbosity_;
170 std::string linsolver = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolver);
171 bdaBridge_ = std::make_unique<detail::BdaSolverInfo<Matrix,Vector>>(accelerator_mode,
172 linear_solver_verbosity,
173 maxit,
174 tolerance,
175 platformID,
176 deviceID,
177 opencl_ilu_parallel,
178 linsolver);
179 }
180
181 void prepare(const Matrix& M, Vector& b)
182 {
183 OPM_TIMEBLOCK(prepare);
184 [[maybe_unused]] const bool firstcall = (this->matrix_ == nullptr);
185 ParentType::prepare(M,b);
186
187#if HAVE_OPENCL
188 // update matrix entries for solvers.
189 if (firstcall && bdaBridge_) {
190 // ebos will not change the matrix object. Hence simply store a pointer
191 // to the original one with a deleter that does nothing.
192 // Outch! We need to be able to scale the linear system! Hence const_cast
193 // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver)
194 bdaBridge_->numJacobiBlocks_ = EWOMS_GET_PARAM(TypeTag, int, NumJacobiBlocks);
195 bdaBridge_->prepare(this->simulator_.vanguard().grid(),
196 this->simulator_.vanguard().cartesianIndexMapper(),
197 this->simulator_.vanguard().schedule().getWellsatEnd(),
198 this->simulator_.vanguard().cellPartition(),
199 this->getMatrix().nonzeroes(), this->useWellConn_);
200 }
201#endif
202 }
203
204
205 void setResidual(Vector& /* b */)
206 {
207 // rhs_ = &b; // Must be handled in prepare() instead.
208 }
209
210 void getResidual(Vector& b) const
211 {
212 b = *(this->rhs_);
213 }
214
215 void setMatrix(const SparseMatrixAdapter& /* M */)
216 {
217 // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
218 }
219
220 bool solve(Vector& x)
221 {
222 if (!bdaBridge_) {
223 return ParentType::solve(x);
224 }
225
226 OPM_TIMEBLOCK(istlSolverEbosBdaSolve);
227 this->solveCount_ += 1;
228 // Write linear system if asked for.
229 const int verbosity = this->prm_[this->activeSolverNum_].template get<int>("verbosity", 0);
230 const bool write_matrix = verbosity > 10;
231 if (write_matrix) {
232 Helper::writeSystem(this->simulator_, //simulator is only used to get names
233 this->getMatrix(),
234 *(this->rhs_),
235 this->comm_.get());
236 }
237
238 // Solve system.
239 Dune::InverseOperatorResult result;
240
241 std::function<void(WellContributions&)> getContribs =
242 [this](WellContributions& w)
243 {
244 this->simulator_.problem().wellModel().getWellContributions(w);
245 };
246 if (!bdaBridge_->apply(*(this->rhs_), this->useWellConn_, getContribs,
247 this->simulator_.gridView().comm().rank(),
248 const_cast<Matrix&>(this->getMatrix()),
249 x, result))
250 {
251 if(bdaBridge_->gpuActive()){
252 // bda solve fails use istl solver setup need to be done since it is not setup in prepare
253 ParentType::prepareFlexibleSolver();
254 }
255 assert(this->flexibleSolver_[this->activeSolverNum_].solver_);
256 this->flexibleSolver_[this->activeSolverNum_].solver_->apply(x, *(this->rhs_), result);
257 }
258
259 // Check convergence, iterations etc.
260 this->checkConvergence(result);
261
262 return this->converged_;
263 }
264
265protected:
266 std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge_;
267}; // end ISTLSolver
268
269} // namespace Opm
270
271#endif // OPM_ISTLSOLVER_EBOS_BDA_HEADER_INCLUDED
Definition findOverlapRowsAndColumns.hpp:31
Interface class adding the update() method to the preconditioner interface.
Definition PreconditionerWithUpdate.hpp:32
Definition AquiferInterface.hpp:35
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition ISTLSolverEbosBda.hpp:98
ISTLSolverEbosBda(const Simulator &simulator)
Construct a system solver.
Definition ISTLSolverEbosBda.hpp:139
ISTLSolverEbosBda(const Simulator &simulator, const FlowLinearSolverParameters &parameters)
Construct a system solver.
Definition ISTLSolverEbosBda.hpp:131
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition ISTLSolverEbos.hpp:144
This class serves to eliminate the need to include the WellContributions into the matrix (with –matri...
Definition WellContributions.hpp:52
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 ISTLSolverEbosBda.hpp:39