StandardPreconditioners_serial.hpp
Go to the documentation of this file.
1/*
2 Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
3 Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
4
5 This file is part of the Open Porous Media project (OPM).
6
7 OPM is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11
12 OPM is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with OPM. If not, see <http://www.gnu.org/licenses/>.
19*/
20
21#ifndef OPM_STANDARDPRECONDITIONERS_SERIAL_HPP
22#define OPM_STANDARDPRECONDITIONERS_SERIAL_HPP
23
24#if HAVE_CUDA
25#include <opm/simulators/linalg/gpuistl/PreconditionerCPUMatrixToGPUMatrix.hpp>
26#endif
27
28namespace Opm {
29
30
31template <class Operator>
32struct StandardPreconditioners<Operator, Dune::Amg::SequentialInformation, typename std::enable_if_t<!Opm::is_gpu_operator_v<Operator>>>
33{
34 static void add()
35 {
36 using namespace Dune;
37 using O = Operator;
38 using C = Dune::Amg::SequentialInformation;
39 using F = PreconditionerFactory<O, C>;
40 using M = typename F::Matrix;
41 using V = typename F::Vector;
42 using P = PropertyTree;
43 F::addCreator("ILU0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
44 const double w = prm.get<double>("relaxation", 1.0);
45 return std::make_shared<ParallelOverlappingILU0<M, V, V, C>>(
46 op.getmat(), 0, w, MILU_VARIANT::ILU);
47 });
48 F::addCreator("DuneILU", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
49 const double w = prm.get<double>("relaxation", 1.0);
50 const int n = prm.get<int>("ilulevel", 0);
51 const bool resort = prm.get<bool>("resort", false);
52 return getRebuildOnUpdateWrapper<Dune::SeqILU<M, V, V>>(op.getmat(), n, w, resort);
53 });
54 F::addCreator("ParOverILU0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
55 const double w = prm.get<double>("relaxation", 1.0);
56 const int n = prm.get<int>("ilulevel", 0);
57 return std::make_shared<ParallelOverlappingILU0<M, V, V, C>>(
58 op.getmat(), n, w, MILU_VARIANT::ILU);
59 });
60 F::addCreator("ILUn", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
61 const int n = prm.get<int>("ilulevel", 0);
62 const double w = prm.get<double>("relaxation", 1.0);
63 return std::make_shared<ParallelOverlappingILU0<M, V, V, C>>(
64 op.getmat(), n, w, MILU_VARIANT::ILU);
65 });
66 F::addCreator("DILU", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
67 DUNE_UNUSED_PARAMETER(prm);
68 return std::make_shared<MultithreadDILU<M, V, V>>(op.getmat());
69 });
70 F::addCreator("Jac", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
71 const int n = prm.get<int>("repeats", 1);
72 const double w = prm.get<double>("relaxation", 1.0);
73 return getDummyUpdateWrapper<SeqJac<M, V, V>>(op.getmat(), n, w);
74 });
75 F::addCreator("GS", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
76 const int n = prm.get<int>("repeats", 1);
77 const double w = prm.get<double>("relaxation", 1.0);
78 return getDummyUpdateWrapper<SeqGS<M, V, V>>(op.getmat(), n, w);
79 });
80 F::addCreator("SOR", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
81 const int n = prm.get<int>("repeats", 1);
82 const double w = prm.get<double>("relaxation", 1.0);
83 return getDummyUpdateWrapper<SeqSOR<M, V, V>>(op.getmat(), n, w);
84 });
85 F::addCreator("SSOR", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
86 const int n = prm.get<int>("repeats", 1);
87 const double w = prm.get<double>("relaxation", 1.0);
88 return getDummyUpdateWrapper<SeqSSOR<M, V, V>>(op.getmat(), n, w);
89 });
90
91 // Only add AMG preconditioners to the factory if the operator
92 // is an actual matrix operator.
93 if constexpr (std::is_same_v<O, Dune::MatrixAdapter<M, V, V>>) {
94 F::addCreator("amg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
95 const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
96 if (smoother == "ILU0" || smoother == "ParOverILU0") {
97 using Smoother = SeqILU<M, V, V>;
98 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
99 } else if (smoother == "Jac") {
100 using Smoother = SeqJac<M, V, V>;
101 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
102 } else if (smoother == "GS") {
103 using Smoother = SeqGS<M, V, V>;
104 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
105 } else if (smoother == "DILU") {
106 using Smoother = MultithreadDILU<M, V, V>;
107 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
108 } else if (smoother == "SOR") {
109 using Smoother = SeqSOR<M, V, V>;
110 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
111 } else if (smoother == "SSOR") {
112 using Smoother = SeqSSOR<M, V, V>;
113 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
114 } else if (smoother == "ILUn") {
115 using Smoother = SeqILU<M, V, V>;
116 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
117 } else {
118 OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + ".");
119 }
120 });
121 F::addCreator("kamg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
122 const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
123 if (smoother == "ILU0" || smoother == "ParOverILU0") {
124 using Smoother = SeqILU<M, V, V>;
125 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
126 } else if (smoother == "Jac") {
127 using Smoother = SeqJac<M, V, V>;
128 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
129 } else if (smoother == "SOR") {
130 using Smoother = SeqSOR<M, V, V>;
131 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
132 } else if (smoother == "GS") {
133 using Smoother = SeqGS<M, V, V>;
134 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
135 } else if (smoother == "SSOR") {
136 using Smoother = SeqSSOR<M, V, V>;
137 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
138 } else if (smoother == "ILUn") {
139 using Smoother = SeqILU<M, V, V>;
140 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
141 } else {
142 OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + ".");
143 }
144 });
145 F::addCreator("famg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
146 if constexpr (std::is_same_v<typename V::field_type, float>) {
147 OPM_THROW(std::logic_error, "famg requires UMFPack which is not available for floats");
148 return nullptr;
149 } else {
150 auto crit = AMGHelper<O, C, M, V>::criterion(prm);
151 Dune::Amg::Parameters parms;
152 parms.setNoPreSmoothSteps(1);
153 parms.setNoPostSmoothSteps(1);
154 return getRebuildOnUpdateWrapper<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
155 }
156 });
157
158#if HAVE_AMGX
159 // Only add AMGX for scalar matrices
160 if constexpr (M::block_type::rows == 1 && M::block_type::cols == 1) {
161 F::addCreator("amgx", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
162 auto prm_copy = prm;
163 prm_copy.put("setup_frequency", Opm::Parameters::Get<Opm::Parameters::CprReuseInterval>());
164 return std::make_shared<Amgx::AmgxPreconditioner<M, V, V>>(op.getmat(), prm_copy);
165 });
166 }
167#endif
168
169#if HAVE_HYPRE
170 // Only add Hypre for scalar matrices
171 if constexpr (M::block_type::rows == 1 && M::block_type::cols == 1 &&
172 std::is_same_v<HYPRE_Real, typename V::field_type>) {
173 F::addCreator("hypre", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
174 return std::make_shared<Hypre::HyprePreconditioner<M, V, V>>(op.getmat(), prm);
175 });
176 }
177#endif
178 }
179
180 // Add CPRW only for the WellModelMatrixAdapter, as the method requires that the operator
181 // has the addWellPressureEquations() method (and a few more) it can not be combined with
182 // a well-less operator such as Dune::MatrixAdapter. For OPM Flow this corresponds to
183 // requiring --matrix-add-well-contributions=false (which is the default).
184 if constexpr (std::is_same_v<O, WellModelMatrixAdapter<M, V, V>>) {
185 F::addCreator(
186 "cprw",
187 [](const O& op, const P& prm, const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
188 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
189 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
190 }
191 using Scalar = typename V::field_type;
192 using LevelTransferPolicy
193 = PressureBhpTransferPolicy<O, Dune::Amg::SequentialInformation, Scalar, false>;
194 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
195 op, prm, weightsCalculator, pressureIndex);
196 });
197 }
198
199 F::addCreator(
200 "cpr",
201 [](const O& op, const P& prm, const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
202 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
203 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
204 }
205 using Scalar = typename V::field_type;
206 using LevelTransferPolicy = PressureTransferPolicy<O, Dune::Amg::SequentialInformation, Scalar, false>;
207 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
208 op, prm, weightsCalculator, pressureIndex);
209 });
210 F::addCreator(
211 "cprt",
212 [](const O& op, const P& prm, const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
213 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
214 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
215 }
216 using Scalar = typename V::field_type;
217 using LevelTransferPolicy = PressureTransferPolicy<O, Dune::Amg::SequentialInformation, Scalar, true>;
218 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
219 op, prm, weightsCalculator, pressureIndex);
220 });
221
222#if HAVE_CUDA
223 // Here we create the *wrapped* GPU preconditioners
224 // meaning they will act as CPU preconditioners on the outside,
225 // but copy data back and forth to the GPU as needed.
226
227 // TODO: Make this use the GPU preconditioner factory once that is up and running.
228 F::addCreator("GPUILU0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
229 const double w = prm.get<double>("relaxation", 1.0);
230 using field_type = typename V::field_type;
231 using GpuILU0 = typename gpuistl::
232 GpuSeqILU0<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
233 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, GpuILU0>>(
234 std::make_shared<GpuILU0>(op.getmat(), w));
235 });
236
237 F::addCreator("GPUILU0Float", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
238 const double w = prm.get<double>("relaxation", 1.0);
239 using block_type = typename V::block_type;
240 using VTo = Dune::BlockVector<Dune::FieldVector<float, block_type::dimension>>;
241 using matrix_type_to =
242 typename Dune::BCRSMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
243 using GpuILU0 = typename gpuistl::
244 GpuSeqILU0<matrix_type_to, gpuistl::GpuVector<float>, gpuistl::GpuVector<float>>;
245 using Adapter = typename gpuistl::PreconditionerAdapter<VTo, VTo, GpuILU0>;
246 using Converter = typename gpuistl::PreconditionerConvertFieldTypeAdapter<Adapter, M, V, V>;
247 auto converted = std::make_shared<Converter>(op.getmat());
248 auto adapted = std::make_shared<Adapter>(std::make_shared<GpuILU0>(converted->getConvertedMatrix(), w));
249 converted->setUnderlyingPreconditioner(adapted);
250 return converted;
251 });
252
253 F::addCreator("GPUJAC", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
254 const double w = prm.get<double>("relaxation", 1.0);
255 using field_type = typename V::field_type;
256 using GPUJac =
257 typename gpuistl::GpuJac<gpuistl::GpuSparseMatrix<field_type>,
258 gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
259
260 using MatrixOwner = Opm::gpuistl::PreconditionerCPUMatrixToGPUMatrix<gpuistl::GpuVector<field_type>,
261 gpuistl::GpuVector<field_type>, GPUJac, M>;
262 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, MatrixOwner>>(
263 std::make_shared<MatrixOwner>(op.getmat(), w));
264 });
265
266 F::addCreator("OPMGPUILU0", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t) {
267 const bool split_matrix = prm.get<bool>("split_matrix", true);
268 const bool tune_gpu_kernels = prm.get<bool>("tune_gpu_kernels", true);
269 const int mixed_precision_scheme = prm.get<int>("mixed_precision_scheme", 0);
270
271 using field_type = typename V::field_type;
272 using GPUILU0 = typename gpuistl::OpmGpuILU0<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
273 using MatrixOwner = Opm::gpuistl::PreconditionerCPUMatrixToGPUMatrix<gpuistl::GpuVector<field_type>,
274 gpuistl::GpuVector<field_type>, GPUILU0, M>;
275 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, MatrixOwner>>(
276 // Note: op.getmat() is passed twice, because the ILU0 needs both the CPU and GPU matrix.
277 // The first argument will be converted to a GPU matrix, and the second one is used as a CPU matrix.
278 std::make_shared<MatrixOwner>(op.getmat(), op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme));
279 });
280
281
282 F::addCreator("GPUDILU", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t) {
283 const bool split_matrix = prm.get<bool>("split_matrix", true);
284 const bool tune_gpu_kernels = prm.get<bool>("tune_gpu_kernels", true);
285 const int mixed_precision_scheme = prm.get<int>("mixed_precision_scheme", 0);
286 using field_type = typename V::field_type;
287 using GPUDILU = typename gpuistl::GpuDILU<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
288 using MatrixOwner = Opm::gpuistl::PreconditionerCPUMatrixToGPUMatrix<gpuistl::GpuVector<field_type>,
289 gpuistl::GpuVector<field_type>, GPUDILU, M>;
290 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, MatrixOwner>>(
291 // Note: op.getmat() is passed twice, because the DILU needs both the CPU and GPU matrix.
292 // The first argument will be converted to a GPU matrix, and the second one is used as a CPU matrix.
293 std::make_shared<MatrixOwner>(op.getmat(), op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme));
294 });
295
296 F::addCreator("GPUDILUFloat", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t) {
297 const bool split_matrix = prm.get<bool>("split_matrix", true);
298 const bool tune_gpu_kernels = prm.get<bool>("tune_gpu_kernels", true);
299 const int mixed_precision_scheme = prm.get<int>("mixed_precision_scheme", 0);
300
301 using block_type = typename V::block_type;
302 using VTo = Dune::BlockVector<Dune::FieldVector<float, block_type::dimension>>;
303 using matrix_type_to = typename Dune::BCRSMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
304 using GpuDILU = typename gpuistl::GpuDILU<matrix_type_to, gpuistl::GpuVector<float>, gpuistl::GpuVector<float>>;
305 using MatrixOwner = Opm::gpuistl::PreconditionerCPUMatrixToGPUMatrix<gpuistl::GpuVector<float>,
306 gpuistl::GpuVector<float>, GpuDILU, matrix_type_to>;
307 using Adapter = typename gpuistl::PreconditionerAdapter<VTo, VTo, MatrixOwner>;
308 using Converter = typename gpuistl::PreconditionerConvertFieldTypeAdapter<Adapter, M, V, V>;
309
310
311 auto converted = std::make_shared<Converter>(op.getmat());
312 // Note: converted->getConvertedMatrix() is passed twice, because the DILU needs both the CPU and GPU matrix.
313 // The first argument will be converted to a GPU matrix, and the second one is used as a CPU matrix.
314 auto adapted = std::make_shared<Adapter>(std::make_shared<MatrixOwner>(
315 converted->getConvertedMatrix(), converted->getConvertedMatrix(),
316 split_matrix, tune_gpu_kernels, mixed_precision_scheme));
317 converted->setUnderlyingPreconditioner(adapted);
318 return converted;
319 });
320#endif // HAVE_CUDA
321 }
322};
323
324
325} // namespace Opm
326
327
328#endif // OPM_STANDARDPRECONDITIONERS_SERIAL_HPP
Dune::MultithreadDILU
The OpenMP thread parallelized DILU preconditioner.
Definition: DILU.hpp:53
Opm::PreconditionerFactory
Definition: PreconditionerFactory.hpp:64
Opm::PressureBhpTransferPolicy
Definition: PressureBhpTransferPolicy.hpp:99
Opm::PressureTransferPolicy
Definition: PressureTransferPolicy.hpp:55
Opm::PropertyTree
Hierarchical collection of key/value pairs.
Definition: PropertyTree.hpp:39
Opm::gpuistl::GpuDILU
DILU preconditioner on the GPU.
Definition: GpuDILU.hpp:50
Opm::gpuistl::GpuJac
Jacobi preconditioner on the GPU.
Definition: GpuJac.hpp:47
Opm::gpuistl::OpmGpuILU0
ILU0 preconditioner on the GPU.
Definition: OpmGpuILU0.hpp:51
Opm::gpuistl::PreconditionerAdapter
Makes a CUDA preconditioner available to a CPU simulator.
Definition: PreconditionerAdapter.hpp:43
Opm::gpuistl::PreconditionerCPUMatrixToGPUMatrix
Convert a CPU matrix to a GPU matrix and use a CUDA preconditioner on the GPU.
Definition: PreconditionerCPUMatrixToGPUMatrix.hpp:42
Opm::gpuistl::PreconditionerConvertFieldTypeAdapter
Converts the field type (eg. double to float) to benchmark single precision preconditioners.
Definition: PreconditionerConvertFieldTypeAdapter.hpp:86
Dune
Definition: fvbaseprimaryvariables.hh:141
Opm
Definition: blackoilboundaryratevector.hh:39
Opm::MILU_VARIANT::ILU
@ ILU
Do not perform modified ILU.
Opm::AMGHelper
Definition: PreconditionerFactory.hpp:43
Opm::AMGHelper::criterion
static Criterion criterion(const PropertyTree &prm)
Definition: StandardPreconditioners_mpi.hpp:83
Opm::StandardPreconditioners
Definition: StandardPreconditioners_mpi.hpp:128