My Project
Loading...
Searching...
No Matches
PreconditionerConvertFieldTypeAdapter.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_PRECONDITIONERCONVERTOFLOATADAPTER_HPP
20#define OPM_PRECONDITIONERCONVERTOFLOATADAPTER_HPP
21#include <cusparse.h>
22#include <dune/istl/bcrsmatrix.hh>
23#include <dune/istl/preconditioner.hh>
24#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
25#include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
26#include <opm/simulators/linalg/gpuistl/detail/CuMatrixDescription.hpp>
27#include <opm/simulators/linalg/gpuistl/detail/CuSparseHandle.hpp>
28#include <opm/simulators/linalg/gpuistl/detail/CuSparseResource.hpp>
29#include <opm/simulators/linalg/gpuistl/detail/cusparse_constants.hpp>
30#include <opm/simulators/linalg/gpuistl/detail/cusparse_safe_call.hpp>
31#include <opm/simulators/linalg/gpuistl/detail/preconditioner_should_call_post_pre.hpp>
32
33
34namespace Opm::gpuistl
35{
84template <class CudaPreconditionerType, class M, class X, class Y, int l = 1>
86{
87public:
89 using matrix_type = typename std::remove_const<M>::type;
90
92 using domain_type = X;
94 using range_type = Y;
96 using field_type = typename X::field_type;
97
98
99 using domain_type_to = typename CudaPreconditionerType::domain_type;
101 using range_type_to = typename CudaPreconditionerType::range_type;
103 using field_type_to = typename domain_type_to::field_type;
104
105
106 using block_type = typename domain_type::block_type;
107
108 using XTo = Dune::BlockVector<Dune::FieldVector<field_type_to, block_type::dimension>>;
109 using YTo = Dune::BlockVector<Dune::FieldVector<field_type_to, block_type::dimension>>;
110 using matrix_type_to =
111 typename Dune::BCRSMatrix<Dune::FieldMatrix<field_type_to, block_type::dimension, block_type::dimension>>;
112
121 : m_matrix(matrix)
122 , m_convertedMatrix(createConvertedMatrix())
123 {
124 }
125
127 virtual void pre([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
128 {
130 "We currently do not support Preconditioner::pre().");
131 }
132
134 virtual void apply(X& v, const Y& d) override
135 {
136 OPM_ERROR_IF(!m_underlyingPreconditioner,
137 "You need to set the underlying preconditioner with setUnderlyingPreconditioner.");
138 XTo convertedV(v.N());
139 for (size_t i = 0; i < v.N(); ++i) {
140 for (size_t j = 0; j < block_type::dimension; ++j) {
141 // This is probably unnecessary, but doing it anyway:
142 convertedV[i][j] = field_type_to(v[i][j]);
143 }
144 }
145 YTo convertedD(d.N());
146 for (size_t i = 0; i < d.N(); ++i) {
147 for (size_t j = 0; j < block_type::dimension; ++j) {
148 convertedD[i][j] = field_type_to(d[i][j]);
149 }
150 }
151
152 m_underlyingPreconditioner->apply(convertedV, convertedD);
153
154 for (size_t i = 0; i < v.N(); ++i) {
155 for (size_t j = 0; j < block_type::dimension; ++j) {
156 v[i][j] = field_type(convertedV[i][j]);
157 }
158 }
159 }
160
162 virtual void post([[maybe_unused]] X& x) override
163 {
165 "We currently do not support Preconditioner::post().");
166 }
167
169 virtual Dune::SolverCategory::Category category() const override
170 {
171 return m_underlyingPreconditioner->category();
172 }
173
174 virtual void update() override
175 {
176 OPM_ERROR_IF(!m_underlyingPreconditioner,
177 "You need to set the underlying preconditioner with setUnderlyingPreconditioner.");
178 updateMatrix();
179 m_underlyingPreconditioner->update();
180 }
181
182 const matrix_type_to& getConvertedMatrix() const
183 {
184 return m_convertedMatrix;
185 }
186
187 void setUnderlyingPreconditioner(const std::shared_ptr<CudaPreconditionerType>& conditioner)
188 {
189 m_underlyingPreconditioner = conditioner;
190 }
191
192 virtual bool hasPerfectUpdate() const override {
193 return m_underlyingPreconditioner->hasPerfectUpdate();
194 }
195
196
197private:
198 void updateMatrix()
199 {
200 const auto nnz = m_matrix.nonzeroes() * m_matrix[0][0].N() * m_matrix[0][0].N();
201 const auto dataPointerIn = static_cast<const field_type*>(&((m_matrix[0][0][0][0])));
202 auto dataPointerOut = static_cast<field_type_to*>(&((m_convertedMatrix[0][0][0][0])));
203
204 std::vector<field_type_to> buffer(nnz, 0);
205 for (size_t i = 0; i < nnz; ++i) {
206 dataPointerOut[i] = field_type_to(dataPointerIn[i]);
207 }
208 }
209 matrix_type_to createConvertedMatrix()
210 {
211 // TODO: Check if this whole conversion can be done more efficiently.
212 const auto N = m_matrix.N();
213 matrix_type_to matrixBuilder(N, N, m_matrix.nonzeroes(), matrix_type_to::row_wise);
214 {
215 auto rowIn = m_matrix.begin();
216 for (auto rowOut = matrixBuilder.createbegin(); rowOut != matrixBuilder.createend(); ++rowOut) {
217 for (auto column = rowIn->begin(); column != rowIn->end(); ++column) {
218 rowOut.insert(column.index());
219 }
220 ++rowIn;
221 }
222 }
223
224 for (auto row = m_matrix.begin(); row != m_matrix.end(); ++row) {
225 for (auto column = row->begin(); column != row->end(); ++column) {
226 for (size_t i = 0; i < block_type::dimension; ++i) {
227 for (size_t j = 0; j < block_type::dimension; ++j) {
228 matrixBuilder[row.index()][column.index()][i][j]
229 = field_type_to(m_matrix[row.index()][column.index()][i][j]);
230 }
231 }
232 }
233 }
234
235 return matrixBuilder;
236 }
237 const M& m_matrix;
238 matrix_type_to m_convertedMatrix;
240 std::shared_ptr<CudaPreconditionerType> m_underlyingPreconditioner;
241};
242} // end namespace Opm::gpuistl
243
244#endif
Interface class adding the update() method to the preconditioner interface.
Definition PreconditionerWithUpdate.hpp:32
Converts the field type (eg.
Definition PreconditionerConvertFieldTypeAdapter.hpp:86
virtual void post(X &x) override
Not used at the moment.
Definition PreconditionerConvertFieldTypeAdapter.hpp:162
typename X::field_type field_type
The field type of the preconditioner.
Definition PreconditionerConvertFieldTypeAdapter.hpp:96
X domain_type
The domain type of the preconditioner.
Definition PreconditionerConvertFieldTypeAdapter.hpp:92
virtual void pre(X &x, Y &b) override
Not used at the moment.
Definition PreconditionerConvertFieldTypeAdapter.hpp:127
PreconditionerConvertFieldTypeAdapter(const M &matrix)
Constructor.
Definition PreconditionerConvertFieldTypeAdapter.hpp:120
virtual void apply(X &v, const Y &d) override
Apply the preconditoner.
Definition PreconditionerConvertFieldTypeAdapter.hpp:134
typename domain_type_to::field_type field_type_to
The field type of the preconditioner.
Definition PreconditionerConvertFieldTypeAdapter.hpp:103
typename std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition PreconditionerConvertFieldTypeAdapter.hpp:89
typename CudaPreconditionerType::range_type range_type_to
The range type of the preconditioner.
Definition PreconditionerConvertFieldTypeAdapter.hpp:101
virtual Dune::SolverCategory::Category category() const override
Category of the preconditioner (see SolverCategory::Category)
Definition PreconditionerConvertFieldTypeAdapter.hpp:169
Y range_type
The range type of the preconditioner.
Definition PreconditionerConvertFieldTypeAdapter.hpp:94
constexpr bool shouldCallPreconditionerPost()
Tests (compile time) if the preconditioner type needs to call post() after a call to apply(....
Definition preconditioner_should_call_post_pre.hpp:52
constexpr bool shouldCallPreconditionerPre()
Tests (compile time) if the preconditioner type needs to call pre() before a call to apply()
Definition preconditioner_should_call_post_pre.hpp:34