My Project
Loading...
Searching...
No Matches
SolverAdapter.hpp
1/*
2 Copyright 2022-2023 SINTEF AS
3
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#ifndef OPM_SOLVERADAPTER_HPP
20#define OPM_SOLVERADAPTER_HPP
21
22#include <memory>
23
24#include <dune/istl/bcrsmatrix.hh>
25#include <dune/istl/operators.hh>
26#include <dune/istl/schwarz.hh>
27#include <dune/istl/solver.hh>
28#include <opm/common/ErrorMacros.hpp>
29#include <opm/simulators/linalg/cuistl/CuBlockPreconditioner.hpp>
30#include <opm/simulators/linalg/cuistl/CuOwnerOverlapCopy.hpp>
31#include <opm/simulators/linalg/cuistl/CuSparseMatrix.hpp>
32#include <opm/simulators/linalg/cuistl/CuVector.hpp>
33#include <opm/simulators/linalg/cuistl/PreconditionerAdapter.hpp>
34#include <opm/simulators/linalg/cuistl/detail/has_function.hpp>
35
36
37
38namespace Opm::cuistl
39{
45template <class Operator, template <class> class UnderlyingSolver, class X>
46class SolverAdapter : public Dune::IterativeSolver<X, X>
47{
48public:
49 using typename Dune::IterativeSolver<X, X>::domain_type;
50 using typename Dune::IterativeSolver<X, X>::range_type;
51 using typename Dune::IterativeSolver<X, X>::field_type;
52 using typename Dune::IterativeSolver<X, X>::real_type;
53 using typename Dune::IterativeSolver<X, X>::scalar_real_type;
54 static constexpr auto block_size = domain_type::block_type::dimension;
55 using XGPU = Opm::cuistl::CuVector<real_type>;
56
57 // TODO: Use a std::forward
58 SolverAdapter(Operator& op,
59 Dune::ScalarProduct<X>& sp,
60 std::shared_ptr<Dune::Preconditioner<X, X>> prec,
61 scalar_real_type reduction,
62 int maxit,
63 int verbose)
64 : Dune::IterativeSolver<X, X>(op, sp, *prec, reduction, maxit, verbose)
65 , m_opOnCPUWithMatrix(op)
66 , m_matrix(CuSparseMatrix<real_type>::fromMatrix(op.getmat()))
67 , m_underlyingSolver(constructSolver(prec, reduction, maxit, verbose))
68 {
69 }
70
71 virtual void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override
72 {
73 // TODO: Can we do this without reimplementing the other function?
74 // TODO: [perf] Do we need to update the matrix every time? Probably yes
75 m_matrix.updateNonzeroValues(m_opOnCPUWithMatrix.getmat());
76
77 if (!m_inputBuffer) {
78 m_inputBuffer.reset(new XGPU(b.dim()));
79 m_outputBuffer.reset(new XGPU(x.dim()));
80 }
81
82 m_inputBuffer->copyFromHost(b);
83 // TODO: [perf] do we need to copy x here?
84 m_outputBuffer->copyFromHost(x);
85
86 m_underlyingSolver.apply(*m_outputBuffer, *m_inputBuffer, reduction, res);
87
88 // TODO: [perf] do we need to copy b here?
89 m_inputBuffer->copyToHost(b);
90 m_outputBuffer->copyToHost(x);
91 }
92 virtual void apply(X& x, X& b, Dune::InverseOperatorResult& res) override
93 {
94 // TODO: [perf] Do we need to update the matrix every time? Probably yes
95 m_matrix.updateNonzeroValues(m_opOnCPUWithMatrix.getmat());
96
97 if (!m_inputBuffer) {
98 m_inputBuffer.reset(new XGPU(b.dim()));
99 m_outputBuffer.reset(new XGPU(x.dim()));
100 }
101
102 m_inputBuffer->copyFromHost(b);
103 // TODO: [perf] do we need to copy x here?
104 m_outputBuffer->copyFromHost(x);
105
106 m_underlyingSolver.apply(*m_outputBuffer, *m_inputBuffer, res);
107
108 // TODO: [perf] do we need to copy b here?
109 m_inputBuffer->copyToHost(b);
110 m_outputBuffer->copyToHost(x);
111 }
112
113private:
114 Operator& m_opOnCPUWithMatrix;
116
117 UnderlyingSolver<XGPU> m_underlyingSolver;
118
119
120 // TODO: Use a std::forward
121 UnderlyingSolver<XGPU> constructSolver(std::shared_ptr<Dune::Preconditioner<X, X>> prec,
122 scalar_real_type reduction,
123 int maxit,
124 int verbose)
125 {
126
129 "Currently we only support operators of type MatrixAdapter in the CUDA solver. "
130 "Use --matrix-add-well-contributions=true. "
131 "Using WellModelMatrixAdapter with SolverAdapter is not well-defined so it will not work well, or at all.");
132
133
134
135
136
138 // TODO: See the below TODO over the definition of precHolder in the other branch
139 // TODO: We are currently double wrapping preconditioners in the preconditioner factory to be extra
140 // compatible
141 // with CPU. Probably a cleaner way eventually would be to do more modifications to the flexible
142 // solver to accomodate the pure GPU better.
143 auto precAsHolder = std::dynamic_pointer_cast<PreconditionerHolder<X, X>>(prec);
144 if (!precAsHolder) {
145 OPM_THROW(std::invalid_argument,
146 "The preconditioner needs to be a CUDA preconditioner (eg. CuILU0) wrapped in a "
147 "Opm::cuistl::PreconditionerAdapter wrapped in a "
148 "Opm::cuistl::CuBlockPreconditioner. If you are unsure what this means, set "
149 "preconditioner to 'CUILU0'"); // TODO: Suggest a better preconditioner
150 }
151
152 auto preconditionerAdapter = precAsHolder->getUnderlyingPreconditioner();
154 = std::dynamic_pointer_cast<PreconditionerHolder<XGPU, XGPU>>(preconditionerAdapter);
156 OPM_THROW(std::invalid_argument,
157 "The preconditioner needs to be a CUDA preconditioner (eg. CuILU0) wrapped in a "
158 "Opm::cuistl::PreconditionerAdapter wrapped in a "
159 "Opm::cuistl::CuBlockPreconditioner. If you are unsure what this means, set "
160 "preconditioner to 'CUILU0'"); // TODO: Suggest a better preconditioner
161 }
162 // We need to get the underlying preconditioner:
163 auto preconditionerReallyOnGPU = preconditionerAdapterAsHolder->getUnderlyingPreconditioner();
164 const auto& communication = m_opOnCPUWithMatrix.getCommunication();
165
167 using SchwarzOperator
168 = Dune::OverlappingSchwarzOperator<CuSparseMatrix<real_type>, XGPU, XGPU, CudaCommunication>;
169 auto cudaCommunication = std::make_shared<CudaCommunication>(communication);
170
171 auto mpiPreconditioner = std::make_shared<CuBlockPreconditioner<XGPU, XGPU, CudaCommunication>>(
173
174 auto scalarProduct = std::make_shared<Dune::ParallelScalarProduct<XGPU, CudaCommunication>>(
175 cudaCommunication, m_opOnCPUWithMatrix.category());
176
177
178 // NOTE: Ownsership of cudaCommunication is handled by mpiPreconditioner. However, just to make sure we
179 // remember
180 // this, we add this check. So remember that we hold one count in this scope, and one in the
181 // CuBlockPreconditioner instance. We accomedate for the fact that it could be passed around in
182 // CuBlockPreconditioner, hence we do not test for != 2
183 OPM_ERROR_IF(cudaCommunication.use_count() < 2, "Internal error. Shared pointer not owned properly.");
184 auto overlappingCudaOperator = std::make_shared<SchwarzOperator>(m_matrix, *cudaCommunication);
185
188 } else {
189 // TODO: Fix the reliance on casting here. This is a code smell to a certain degree, and assumes
190 // a certain setup beforehand. The only reason we do it this way is to minimize edits to the
191 // flexible solver. We could design it differently, but keep this for the time being until
192 // we figure out how we want to GPU-ify the rest of the system.
193 auto precAsHolder = std::dynamic_pointer_cast<PreconditionerHolder<XGPU, XGPU>>(prec);
194 if (!precAsHolder) {
195 OPM_THROW(std::invalid_argument,
196 "The preconditioner needs to be a CUDA preconditioner wrapped in a "
197 "Opm::cuistl::PreconditionerHolder (eg. CuILU0).");
198 }
199 auto preconditionerOnGPU = precAsHolder->getUnderlyingPreconditioner();
200
201 auto matrixOperator
202 = std::make_shared<Dune::MatrixAdapter<CuSparseMatrix<real_type>, XGPU, XGPU>>(m_matrix);
203 auto scalarProduct = std::make_shared<Dune::SeqScalarProduct<XGPU>>();
206 }
207 }
208
209 std::unique_ptr<XGPU> m_inputBuffer;
210 std::unique_ptr<XGPU> m_outputBuffer;
211};
212} // namespace Opm::cuistl
213
214#endif
Definition AquiferInterface.hpp:35
The CuSparseMatrix class simple wrapper class for a CuSparse matrix.
Definition CuSparseMatrix.hpp:48
void updateNonzeroValues(const MatrixType &matrix, bool copyNonZeroElementsDirectly=false)
updateNonzeroValues updates the non-zero values by using the non-zero values of the supplied matrix
Definition CuSparseMatrix.cpp:139
Wraps a CUDA solver to work with CPU data.
Definition SolverAdapter.hpp:47
The has_communication class checks if the type has the member function getCommunication.
Definition has_function.hpp:96
The is_a_well_operator class tries to guess if the operator is a well type operator.
Definition has_function.hpp:114