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#if USE_HIP
26#include <opm/simulators/linalg/gpuistl_hip/PreconditionerCPUMatrixToGPUMatrix.hpp>
27#else
29#endif
30#endif
31
32#include <functional>
33#include <memory>
34#include <type_traits>
35
36namespace Opm {
37
38template <class X, class Y>
40{
41 public:
43 virtual void update() override {};
44 virtual bool hasPerfectUpdate() const override {return true;}
45 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& y) override {};
46 virtual void post ([[maybe_unused]] X& x) override {};
47 virtual void apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y) override {};
48 virtual Dune::SolverCategory::Category category() const override { return Dune::SolverCategory::sequential; };
49};
50
51
52template <class Operator>
53struct StandardPreconditioners<Operator, Dune::Amg::SequentialInformation, typename std::enable_if_t<!Opm::is_gpu_operator_v<Operator>>>
54{
55 static void add()
56 {
57 using namespace Dune;
58 using O = Operator;
59 using C = Dune::Amg::SequentialInformation;
61 using M = typename F::Matrix;
62 using V = typename F::Vector;
63 using P = PropertyTree;
64 F::addCreator("ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
65 const double w = prm.get<double>("relaxation", 1.0);
66 return std::make_shared<ParallelOverlappingILU0<M, V, V, C>>(
67 op.getmat(), 0, w, MILU_VARIANT::ILU);
68 });
69 F::addCreator("duneilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
70 const double w = prm.get<double>("relaxation", 1.0);
71 const int n = prm.get<int>("ilulevel", 0);
72 const bool resort = prm.get<bool>("resort", false);
73 return getRebuildOnUpdateWrapper<Dune::SeqILU<M, V, V>>(std::cref(op.getmat()), n, w, resort);
74 });
75 F::addCreator("paroverilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
76 const double w = prm.get<double>("relaxation", 1.0);
77 const int n = prm.get<int>("ilulevel", 0);
78 return std::make_shared<ParallelOverlappingILU0<M, V, V, C>>(
79 op.getmat(), n, w, MILU_VARIANT::ILU);
80 });
81 F::addCreator("ilun", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
82 const int n = prm.get<int>("ilulevel", 0);
83 const double w = prm.get<double>("relaxation", 1.0);
84 return std::make_shared<ParallelOverlappingILU0<M, V, V, C>>(
85 op.getmat(), n, w, MILU_VARIANT::ILU);
86 });
87 F::addCreator("dilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
88 DUNE_UNUSED_PARAMETER(prm);
89 return std::make_shared<MultithreadDILU<M, V, V>>(op.getmat());
90 });
91 F::addCreator("mixed-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
92 DUNE_UNUSED_PARAMETER(prm);
93 DUNE_UNUSED_PARAMETER(op);
94 return std::make_shared<TrivialPreconditioner<V,V>>();
95 });
96 F::addCreator("mixed-dilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
97 DUNE_UNUSED_PARAMETER(prm);
98 DUNE_UNUSED_PARAMETER(op);
99 return std::make_shared<TrivialPreconditioner<V,V>>();
100 });
101 F::addCreator("jac", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
102 const int n = prm.get<int>("repeats", 1);
103 const double w = prm.get<double>("relaxation", 1.0);
104 return getDummyUpdateWrapper<SeqJac<M, V, V>>(op.getmat(), n, w);
105 });
106 F::addCreator("gs", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
107 const int n = prm.get<int>("repeats", 1);
108 const double w = prm.get<double>("relaxation", 1.0);
109 return getDummyUpdateWrapper<SeqGS<M, V, V>>(op.getmat(), n, w);
110 });
111 F::addCreator("sor", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
112 const int n = prm.get<int>("repeats", 1);
113 const double w = prm.get<double>("relaxation", 1.0);
114 return getDummyUpdateWrapper<SeqSOR<M, V, V>>(op.getmat(), n, w);
115 });
116 F::addCreator("ssor", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
117 const int n = prm.get<int>("repeats", 1);
118 const double w = prm.get<double>("relaxation", 1.0);
119 return getDummyUpdateWrapper<SeqSSOR<M, V, V>>(op.getmat(), n, w);
120 });
121
122 // Only add AMG preconditioners to the factory if the operator
123 // is an actual matrix operator.
124 if constexpr (std::is_same_v<O, Dune::MatrixAdapter<M, V, V>>) {
125 F::addCreator("amg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
126 std::string smoother = prm.get<std::string>("smoother", "paroverilu0");
127 // Make the smoother type lowercase for internal canonical representation
128 std::ranges::transform(smoother, smoother.begin(), ::tolower);
129 if (smoother == "ilu0" || smoother == "paroverilu0") {
130 using Smoother = SeqILU<M, V, V>;
131 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
132 } else if (smoother == "jac") {
133 using Smoother = SeqJac<M, V, V>;
134 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
135 } else if (smoother == "gs") {
136 using Smoother = SeqGS<M, V, V>;
137 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
138 } else if (smoother == "dilu") {
139 using Smoother = MultithreadDILU<M, V, V>;
140 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
141 } else if (smoother == "sor") {
142 using Smoother = SeqSOR<M, V, V>;
143 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
144 } else if (smoother == "ssor") {
145 using Smoother = SeqSSOR<M, V, V>;
146 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
147 } else if (smoother == "ilun") {
148 using Smoother = SeqILU<M, V, V>;
149 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm);
150 } else {
151 OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + ".");
152 }
153 });
154 F::addCreator("kamg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
155 std::string smoother = prm.get<std::string>("smoother", "paroverilu0");
156 // Make the smoother type lowercase for internal canonical representation
157 std::ranges::transform(smoother, smoother.begin(), ::tolower);
158 if (smoother == "ilu0" || smoother == "paroverilu0") {
159 using Smoother = SeqILU<M, V, V>;
160 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
161 } else if (smoother == "jac") {
162 using Smoother = SeqJac<M, V, V>;
163 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
164 } else if (smoother == "sor") {
165 using Smoother = SeqSOR<M, V, V>;
166 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
167 } else if (smoother == "gs") {
168 using Smoother = SeqGS<M, V, V>;
169 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
170 } else if (smoother == "ssor") {
171 using Smoother = SeqSSOR<M, V, V>;
172 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
173 } else if (smoother == "ilun") {
174 using Smoother = SeqILU<M, V, V>;
175 return AMGHelper<O, C, M, V>::template makeAmgPreconditioner<Smoother>(op, prm, true);
176 } else {
177 OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + ".");
178 }
179 });
180 F::addCreator("famg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
181 if constexpr (std::is_same_v<typename V::field_type, float>) {
182 OPM_THROW(std::logic_error, "famg requires UMFPack which is not available for floats");
183 return nullptr;
184 } else {
185 auto crit = AMGHelper<O, C, M, V>::criterion(prm);
186 Dune::Amg::Parameters parms;
187 parms.setNoPreSmoothSteps(1);
188 parms.setNoPostSmoothSteps(1);
189 return getRebuildOnUpdateWrapper<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
190 }
191 });
192
193#if HAVE_AMGX
194 // Only add AMGX for scalar matrices
195 if constexpr (M::block_type::rows == 1 && M::block_type::cols == 1) {
196 F::addCreator("amgx", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
197 auto prm_copy = prm;
198 prm_copy.put("setup_frequency", Opm::Parameters::Get<Opm::Parameters::CprReuseInterval>());
199 return std::make_shared<Amgx::AmgxPreconditioner<M, V, V>>(op.getmat(), prm_copy);
200 });
201 }
202#endif
203
204#if HAVE_HYPRE
205 // Only add Hypre for scalar matrices
206 if constexpr (M::block_type::rows == 1 && M::block_type::cols == 1 &&
207 std::is_same_v<HYPRE_Real, typename V::field_type>) {
208 F::addCreator("hypre", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
209 return std::make_shared<Hypre::HyprePreconditioner<M, V, V, Dune::Amg::SequentialInformation>>(op.getmat(), prm, Dune::Amg::SequentialInformation());
210 });
211 }
212#endif
213 }
214
215 // Add CPRW only for the WellModelMatrixAdapter, as the method requires that the operator
216 // has the addWellPressureEquations() method (and a few more) it can not be combined with
217 // a well-less operator such as Dune::MatrixAdapter. For OPM Flow this corresponds to
218 // requiring --matrix-add-well-contributions=false (which is the default).
219 if constexpr (std::is_same_v<O, WellModelMatrixAdapter<M, V, V>>) {
220 F::addCreator(
221 "cprw",
222 [](const O& op, const P& prm, const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
223 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
224 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
225 }
226 using Scalar = typename V::field_type;
227 using LevelTransferPolicy
229 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
230 op, prm, weightsCalculator, pressureIndex);
231 });
232 }
233
234 F::addCreator(
235 "cpr",
236 [](const O& op, const P& prm, const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
237 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
238 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
239 }
240 using Scalar = typename V::field_type;
242 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
243 op, prm, weightsCalculator, pressureIndex);
244 });
245 F::addCreator(
246 "cprt",
247 [](const O& op, const P& prm, const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
248 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
249 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
250 }
251 using Scalar = typename V::field_type;
253 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
254 op, prm, weightsCalculator, pressureIndex);
255 });
256
257#if HAVE_CUDA
258 // Here we create the *wrapped* GPU preconditioners
259 // meaning they will act as CPU preconditioners on the outside,
260 // but copy data back and forth to the GPU as needed.
261
262 // TODO: Make this use the GPU preconditioner factory once that is up and running.
263 F::addCreator("gpuilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
264 const double w = prm.get<double>("relaxation", 1.0);
265 using field_type = typename V::field_type;
266 using GpuILU0 = typename gpuistl::
267 GpuSeqILU0<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
268 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, GpuILU0>>(
269 std::make_shared<GpuILU0>(op.getmat(), w));
270 });
271
272 F::addCreator("gpuilu0float", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
273 const double w = prm.get<double>("relaxation", 1.0);
274 using block_type = typename V::block_type;
275 using VTo = Dune::BlockVector<Dune::FieldVector<float, block_type::dimension>>;
276 using matrix_type_to =
277 typename Dune::BCRSMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
278 using GpuILU0 = typename gpuistl::
279 GpuSeqILU0<matrix_type_to, gpuistl::GpuVector<float>, gpuistl::GpuVector<float>>;
282 auto converted = std::make_shared<Converter>(op.getmat());
283 auto adapted = std::make_shared<Adapter>(std::make_shared<GpuILU0>(converted->getConvertedMatrix(), w));
284 converted->setUnderlyingPreconditioner(adapted);
285 return converted;
286 });
287
288 F::addCreator("gpujac", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
289 const double w = prm.get<double>("relaxation", 1.0);
290 using field_type = typename V::field_type;
291 using GPUJac =
294
297 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, MatrixOwner>>(
298 std::make_shared<MatrixOwner>(op.getmat(), w));
299 });
300
301 F::addCreator("opmgpuilu0", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t) {
302 const bool split_matrix = prm.get<bool>("split_matrix", true);
303 const bool tune_gpu_kernels = prm.get<bool>("tune_gpu_kernels", true);
304 const int mixed_precision_scheme = prm.get<int>("mixed_precision_scheme", 0);
305
306 using field_type = typename V::field_type;
310 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, MatrixOwner>>(
311 // Note: op.getmat() is passed twice, because the ILU0 needs both the CPU and GPU matrix.
312 // The first argument will be converted to a GPU matrix, and the second one is used as a CPU matrix.
313 std::make_shared<MatrixOwner>(op.getmat(), op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme));
314 });
315
316
317 F::addCreator("gpudilu", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t) {
318 const bool split_matrix = prm.get<bool>("split_matrix", true);
319 const bool tune_gpu_kernels = prm.get<bool>("tune_gpu_kernels", true);
320 const int mixed_precision_scheme = prm.get<int>("mixed_precision_scheme", 0);
321 const bool reorder = prm.get<bool>("reorder", true);
322 using field_type = typename V::field_type;
326 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, MatrixOwner>>(
327 // Note: op.getmat() is passed twice, because the DILU needs both the CPU and GPU matrix.
328 // The first argument will be converted to a GPU matrix, and the second one is used as a CPU matrix.
329 std::make_shared<MatrixOwner>(op.getmat(), op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme, reorder));
330 });
331
332 F::addCreator("gpudilufloat", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t) {
333 const bool split_matrix = prm.get<bool>("split_matrix", true);
334 const bool tune_gpu_kernels = prm.get<bool>("tune_gpu_kernels", true);
335 const int mixed_precision_scheme = prm.get<int>("mixed_precision_scheme", 0);
336 const bool reorder = prm.get<bool>("reorder", true);
337
338 using block_type = typename V::block_type;
339 using VTo = Dune::BlockVector<Dune::FieldVector<float, block_type::dimension>>;
340 using matrix_type_to = typename Dune::BCRSMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
343 gpuistl::GpuVector<float>, GpuDILU, matrix_type_to>;
346
347
348 auto converted = std::make_shared<Converter>(op.getmat());
349 // Note: converted->getConvertedMatrix() is passed twice, because the DILU needs both the CPU and GPU matrix.
350 // The first argument will be converted to a GPU matrix, and the second one is used as a CPU matrix.
351 auto adapted = std::make_shared<Adapter>(std::make_shared<MatrixOwner>(
352 converted->getConvertedMatrix(), converted->getConvertedMatrix(),
353 split_matrix, tune_gpu_kernels, mixed_precision_scheme, reorder));
354 converted->setUnderlyingPreconditioner(adapted);
355 return converted;
356 });
357#endif // HAVE_CUDA
358 }
359};
360
361
362} // namespace Opm
363
364
365#endif // OPM_STANDARDPRECONDITIONERS_SERIAL_HPP
The OpenMP thread parallelized DILU preconditioner.
Definition: DILU.hpp:53
Interface class adding the update() method to the preconditioner interface.
Definition: PreconditionerWithUpdate.hpp:34
Definition: PreconditionerFactory.hpp:64
Definition: PressureBhpTransferPolicy.hpp:99
Definition: PressureTransferPolicy.hpp:55
Hierarchical collection of key/value pairs.
Definition: PropertyTree.hpp:39
Definition: StandardPreconditioners_serial.hpp:40
virtual bool hasPerfectUpdate() const override
Definition: StandardPreconditioners_serial.hpp:44
virtual void pre(X &x, Y &y) override
Definition: StandardPreconditioners_serial.hpp:45
virtual void post(X &x) override
Definition: StandardPreconditioners_serial.hpp:46
virtual void update() override
Definition: StandardPreconditioners_serial.hpp:43
virtual void apply(X &x, const Y &y) override
Definition: StandardPreconditioners_serial.hpp:47
TrivialPreconditioner()
Definition: StandardPreconditioners_serial.hpp:42
virtual Dune::SolverCategory::Category category() const override
Definition: StandardPreconditioners_serial.hpp:48
DILU preconditioner on the GPU.
Definition: GpuDILU.hpp:53
Jacobi preconditioner on the GPU.
Definition: GpuJac.hpp:47
ILU0 preconditioner on the GPU.
Definition: OpmGpuILU0.hpp:51
Makes a CUDA preconditioner available to a CPU simulator.
Definition: PreconditionerAdapter.hpp:43
Convert a CPU matrix to a GPU matrix and use a CUDA preconditioner on the GPU.
Definition: PreconditionerCPUMatrixToGPUMatrix.hpp:42
Converts the field type (eg. double to float) to benchmark single precision preconditioners.
Definition: PreconditionerConvertFieldTypeAdapter.hpp:86
Definition: fvbaseprimaryvariables.hh:161
Definition: blackoilbioeffectsmodules.hh:45
@ ILU
Do not perform modified ILU.
Definition: PreconditionerFactory.hpp:43
static Criterion criterion(const PropertyTree &prm)
Definition: StandardPreconditioners_mpi.hpp:90
Definition: StandardPreconditioners_mpi.hpp:135