My Project
Loading...
Searching...
No Matches
PreconditionerFactory_impl.hpp
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#include <config.h>
22
23#include <opm/common/ErrorMacros.hpp>
24#include <opm/common/TimingMacros.hpp>
25
26#include <opm/simulators/linalg/PreconditionerFactory.hpp>
27
28#include <opm/simulators/linalg/DILU.hpp>
29#include <opm/simulators/linalg/ExtraSmoothers.hpp>
30#include <opm/simulators/linalg/FlexibleSolver.hpp>
31#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
32#include <opm/simulators/linalg/OwningBlockPreconditioner.hpp>
33#include <opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp>
34#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
35#include <opm/simulators/linalg/PressureBhpTransferPolicy.hpp>
36#include <opm/simulators/linalg/PressureTransferPolicy.hpp>
37#include <opm/simulators/linalg/PropertyTree.hpp>
38#include <opm/simulators/linalg/WellOperators.hpp>
40#include <opm/simulators/linalg/ilufirstelement.hh>
41#include <opm/simulators/linalg/matrixblock.hh>
42
43#include <dune/common/unused.hh>
44#include <dune/istl/owneroverlapcopy.hh>
45#include <dune/istl/paamg/amg.hh>
46#include <dune/istl/paamg/fastamg.hh>
47#include <dune/istl/paamg/kamg.hh>
48#include <dune/istl/preconditioners.hh>
49
50// Include all cuistl/GPU preconditioners inside of this headerfile
51#include <opm/simulators/linalg/PreconditionerFactoryGPUIncludeWrapper.hpp>
52
53#if HAVE_HYPRE
54#include <opm/simulators/linalg/HyprePreconditioner.hpp>
55#endif
56
57#if HAVE_AMGX
58#include <opm/simulators/linalg/AmgxPreconditioner.hpp>
59#endif
60
61#include <cassert>
62
63namespace Opm {
64
65template <class Smoother>
67{
68 static auto args(const PropertyTree& prm)
69 {
70 using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
71 SmootherArgs smootherArgs;
72 smootherArgs.iterations = prm.get<int>("iterations", 1);
73 // smootherArgs.overlap=SmootherArgs::vertex;
74 // smootherArgs.overlap=SmootherArgs::none;
75 // smootherArgs.overlap=SmootherArgs::aggregate;
76 smootherArgs.relaxationFactor = prm.get<double>("relaxation", 1.0);
77 return smootherArgs;
78 }
79};
80
81template <class M, class V, class C>
83{
84 static auto args(const PropertyTree& prm)
85 {
87 using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
88 SmootherArgs smootherArgs;
89 smootherArgs.iterations = prm.get<int>("iterations", 1);
90 const int iluwitdh = prm.get<int>("iluwidth", 0);
92 const MILU_VARIANT milu = convertString2Milu(prm.get<std::string>("milutype", std::string("ilu")));
93 smootherArgs.setMilu(milu);
94 // smootherArgs.overlap=SmootherArgs::vertex;
95 // smootherArgs.overlap=SmootherArgs::none;
96 // smootherArgs.overlap=SmootherArgs::aggregate;
97 smootherArgs.relaxationFactor = prm.get<double>("relaxation", 1.0);
98 return smootherArgs;
99 }
100};
101
102// trailing return type with decltype used for detecting existence of setUseFixedOrder member function by overloading the setUseFixedOrder function
103template <typename C>
104auto setUseFixedOrder(C& criterion, bool booleanValue) -> decltype(criterion.setUseFixedOrder(booleanValue))
105{
106 return criterion.setUseFixedOrder(booleanValue); // Set flag to ensure that the matrices in the AMG hierarchy are constructed with deterministic indices.
107}
108template <typename C>
109void setUseFixedOrder(C&, ...)
110{
111 // do nothing, since the function setUseFixedOrder does not exist yet
112}
113
114template <class Operator, class Comm, class Matrix, class Vector>
115typename AMGHelper<Operator, Comm, Matrix, Vector>::Criterion
116AMGHelper<Operator, Comm, Matrix, Vector>::criterion(const PropertyTree& prm)
117{
118 Criterion criterion(15, prm.get<int>("coarsenTarget", 1200));
119 criterion.setDefaultValuesIsotropic(2);
120 criterion.setAlpha(prm.get<double>("alpha", 0.33));
121 criterion.setBeta(prm.get<double>("beta", 1e-5));
122 criterion.setMaxLevel(prm.get<int>("maxlevel", 15));
123 criterion.setSkipIsolated(prm.get<bool>("skip_isolated", false));
124 criterion.setNoPreSmoothSteps(prm.get<int>("pre_smooth", 1));
125 criterion.setNoPostSmoothSteps(prm.get<int>("post_smooth", 1));
126 criterion.setDebugLevel(prm.get<int>("verbosity", 0));
127 // As the default we request to accumulate data to 1 process always as our matrix
128 // graph might be unsymmetric and hence not supported by the PTScotch/ParMetis
129 // calls in DUNE. Accumulating to 1 skips PTScotch/ParMetis
130 criterion.setAccumulate(static_cast<Dune::Amg::AccumulationMode>(prm.get<int>("accumulate", 1)));
131 criterion.setProlongationDampingFactor(prm.get<double>("prolongationdamping", 1.6));
132 criterion.setMaxDistance(prm.get<int>("maxdistance", 2));
133 criterion.setMaxConnectivity(prm.get<int>("maxconnectivity", 15));
134 criterion.setMaxAggregateSize(prm.get<int>("maxaggsize", 6));
135 criterion.setMinAggregateSize(prm.get<int>("minaggsize", 4));
136 setUseFixedOrder(criterion, true); // If possible, set flag to ensure that the matrices in the AMG hierarchy are constructed with deterministic indices.
137 return criterion;
138}
139
140template <class Operator, class Comm, class Matrix, class Vector>
141template <class Smoother>
142typename AMGHelper<Operator, Comm, Matrix, Vector>::PrecPtr
143AMGHelper<Operator, Comm, Matrix, Vector>::makeAmgPreconditioner(const Operator& op,
144 const PropertyTree& prm,
145 bool useKamg)
146{
147 auto crit = criterion(prm);
148 auto sargs = AMGSmootherArgsHelper<Smoother>::args(prm);
149 if (useKamg) {
151 return std::make_shared<Type>(
152 op, crit, sargs, prm.get<std::size_t>("max_krylov", 1), prm.get<double>("min_reduction", 1e-1));
153 } else {
155 return std::make_shared<Type>(op, crit, sargs);
156 }
157}
158
159template <class Operator, class Comm>
161 static void add()
162 {
163 using namespace Dune;
164 using O = Operator;
165 using C = Comm;
167 using M = typename F::Matrix;
168 using V = typename F::Vector;
169 using P = PropertyTree;
170 F::addCreator("ILU0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
171 return createParILU(op, prm, comm, 0);
172 });
173 F::addCreator("ParOverILU0",
174 [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
175 return createParILU(op, prm, comm, prm.get<int>("ilulevel", 0));
176 });
177 F::addCreator("ILUn", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
178 return createParILU(op, prm, comm, prm.get<int>("ilulevel", 0));
179 });
180 F::addCreator("DuneILU", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
181 const int n = prm.get<int>("ilulevel", 0);
182 const double w = prm.get<double>("relaxation", 1.0);
183 const bool resort = prm.get<bool>("resort", false);
185 comm, op.getmat(), n, w, resort);
186 });
187 F::addCreator("DILU", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
190 });
191 F::addCreator("Jac", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
192 const int n = prm.get<int>("repeats", 1);
193 const double w = prm.get<double>("relaxation", 1.0);
195 });
196 F::addCreator("GS", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
197 const int n = prm.get<int>("repeats", 1);
198 const double w = prm.get<double>("relaxation", 1.0);
200 });
201 F::addCreator("SOR", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
202 const int n = prm.get<int>("repeats", 1);
203 const double w = prm.get<double>("relaxation", 1.0);
205 });
206 F::addCreator("SSOR", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
207 const int n = prm.get<int>("repeats", 1);
208 const double w = prm.get<double>("relaxation", 1.0);
210 });
211
212 // Only add AMG preconditioners to the factory if the operator
213 // is the overlapping schwarz operator. This could be extended
214 // later, but at this point no other operators are compatible
215 // with the AMG hierarchy construction.
216 if constexpr (std::is_same_v<O, Dune::OverlappingSchwarzOperator<M, V, V, C>> ||
217 std::is_same_v<O, Opm::GhostLastMatrixAdapter<M, V, V, C>>) {
218 F::addCreator("amg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
219 using PrecPtr = std::shared_ptr<Dune::PreconditionerWithUpdate<V, V>>;
220 const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
221 // TODO: merge this with ILUn, and possibly simplify the factory to only work with ILU?
222 if (smoother == "ILU0" || smoother == "ParOverILU0") {
226 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
227 return prec;
228 } else if (smoother == "DILU") {
230 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
231 using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
232 SmootherArgs sargs;
234 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
235 return prec;
236 } else if (smoother == "Jac") {
238 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
239 using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
240 SmootherArgs sargs;
242 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
243 return prec;
244 } else if (smoother == "GS") {
246 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
247 using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
248 SmootherArgs sargs;
250 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
251 return prec;
252 } else if (smoother == "SOR") {
254 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
255 using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
256 SmootherArgs sargs;
258 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
259 return prec;
260 } else if (smoother == "SSOR") {
262 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
263 using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
264 SmootherArgs sargs;
266 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
267 return prec;
268 } else if (smoother == "ILUn") {
270 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
271 using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
272 SmootherArgs sargs;
274 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
275 return prec;
276 } else {
277 OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + ".");
278 }
279 });
280 }
281
282 F::addCreator("cpr",
283 [](const O& op,
284 const P& prm,
285 const std::function<V()> weightsCalculator,
286 std::size_t pressureIndex,
287 const C& comm) {
289 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
290 OPM_THROW(std::logic_error,
291 "Pressure index out of bounds. It needs to specified for CPR");
292 }
293 using Scalar = typename V::field_type;
295 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
296 op, prm, weightsCalculator, pressureIndex, comm);
297 });
298 F::addCreator("cprt",
299 [](const O& op,
300 const P& prm,
301 const std::function<V()> weightsCalculator,
302 std::size_t pressureIndex,
303 const C& comm) {
305 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
306 OPM_THROW(std::logic_error,
307 "Pressure index out of bounds. It needs to specified for CPR");
308 }
309 using Scalar = typename V::field_type;
311 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
312 op, prm, weightsCalculator, pressureIndex, comm);
313 });
314
315 if constexpr (std::is_same_v<O, WellModelGhostLastMatrixAdapter<M, V, V, true>>) {
316 F::addCreator("cprw",
317 [](const O& op,
318 const P& prm,
319 const std::function<V()> weightsCalculator,
320 std::size_t pressureIndex,
321 const C& comm) {
323 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
324 OPM_THROW(std::logic_error,
325 "Pressure index out of bounds. It needs to specified for CPR");
326 }
327 using Scalar = typename V::field_type;
329 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
330 op, prm, weightsCalculator, pressureIndex, comm);
331 });
332 }
333
334#if HAVE_CUDA
335 F::addCreator("GPUILU0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
336 const double w = prm.get<double>("relaxation", 1.0);
337 using field_type = typename V::field_type;
338 using GpuILU0 = typename gpuistl::
339 GpuSeqILU0<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
340 auto gpuILU0 = std::make_shared<GpuILU0>(op.getmat(), w);
341
342 auto adapted = std::make_shared<gpuistl::PreconditionerAdapter<V, V, GpuILU0>>(gpuILU0);
343 auto wrapped = std::make_shared<gpuistl::GpuBlockPreconditioner<V, V, Comm>>(adapted, comm);
344 return wrapped;
345 });
346
347 F::addCreator("GPUJAC", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
348 const double w = prm.get<double>("relaxation", 1.0);
349 using field_type = typename V::field_type;
350 using GpuJac =
351 typename gpuistl::GpuJac<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
352 auto gpuJac = std::make_shared<GpuJac>(op.getmat(), w);
353
354 auto adapted = std::make_shared<gpuistl::PreconditionerAdapter<V, V, GpuJac>>(gpuJac);
355 auto wrapped = std::make_shared<gpuistl::GpuBlockPreconditioner<V, V, Comm>>(adapted, comm);
356 return wrapped;
357 });
358
359 F::addCreator("GPUDILU", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
360 const bool split_matrix = prm.get<bool>("split_matrix", true);
361 const bool tune_gpu_kernels = prm.get<bool>("tune_gpu_kernels", true);
362 const int mixed_precision_scheme = prm.get<int>("mixed_precision_scheme", 0);
363 using field_type = typename V::field_type;
364 using GpuDILU = typename gpuistl::GpuDILU<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
365 auto gpuDILU = std::make_shared<GpuDILU>(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme);
366
367 auto adapted = std::make_shared<gpuistl::PreconditionerAdapter<V, V, GpuDILU>>(gpuDILU);
368 auto wrapped = std::make_shared<gpuistl::GpuBlockPreconditioner<V, V, Comm>>(adapted, comm);
369 return wrapped;
370 });
371
372 F::addCreator("OPMGPUILU0", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
373 const bool split_matrix = prm.get<bool>("split_matrix", true);
374 const bool tune_gpu_kernels = prm.get<bool>("tune_gpu_kernels", true);
375 const int mixed_precision_scheme = prm.get<int>("mixed_precision_scheme", 0);
376 using field_type = typename V::field_type;
377 using OpmGpuILU0 = typename gpuistl::OpmGpuILU0<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
378 auto gpuilu0 = std::make_shared<OpmGpuILU0>(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme);
379
380 auto adapted = std::make_shared<gpuistl::PreconditionerAdapter<V, V, OpmGpuILU0>>(gpuilu0);
381 auto wrapped = std::make_shared<gpuistl::GpuBlockPreconditioner<V, V, Comm>>(adapted, comm);
382 return wrapped;
383 });
384#endif
385 }
386
387
389 createParILU(const Operator& op, const PropertyTree& prm, const Comm& comm, const int ilulevel)
390 {
392 using M = typename F::Matrix;
393 using V = typename F::Vector;
394
395 const double w = prm.get<double>("relaxation", 1.0);
396 const bool redblack = prm.get<bool>("redblack", false);
397 const bool reorder_spheres = prm.get<bool>("reorder_spheres", false);
398 // Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner.
399 if (ilulevel == 0) {
400 const std::size_t num_interior = interiorIfGhostLast(comm);
401 assert(num_interior <= op.getmat().N());
402 return std::make_shared<ParallelOverlappingILU0<M, V, V, Comm>>(
404 } else {
405 return std::make_shared<ParallelOverlappingILU0<M, V, V, Comm>>(
407 }
408 }
409
414 static std::size_t interiorIfGhostLast(const Comm& comm)
415 {
416 std::size_t interior_count = 0;
417 std::size_t highest_interior_index = 0;
418 const auto& is = comm.indexSet();
419 for (const auto& ind : is) {
420 if (Comm::OwnerSet::contains(ind.local().attribute())) {
422 highest_interior_index = std::max(highest_interior_index, ind.local().local());
423 }
424 }
426 return interior_count;
427 } else {
428 return is.size();
429 }
430 }
431};
432
433template <class Operator>
434struct StandardPreconditioners<Operator, Dune::Amg::SequentialInformation> {
435 static void add()
436 {
437 using namespace Dune;
438 using O = Operator;
439 using C = Dune::Amg::SequentialInformation;
441 using M = typename F::Matrix;
442 using V = typename F::Vector;
443 using P = PropertyTree;
444 F::addCreator("ILU0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
445 const double w = prm.get<double>("relaxation", 1.0);
446 return std::make_shared<ParallelOverlappingILU0<M, V, V, C>>(
447 op.getmat(), 0, w, MILU_VARIANT::ILU);
448 });
449 F::addCreator("DuneILU", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
450 const double w = prm.get<double>("relaxation", 1.0);
451 const int n = prm.get<int>("ilulevel", 0);
452 const bool resort = prm.get<bool>("resort", false);
454 });
455 F::addCreator("ParOverILU0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
456 const double w = prm.get<double>("relaxation", 1.0);
457 const int n = prm.get<int>("ilulevel", 0);
458 return std::make_shared<ParallelOverlappingILU0<M, V, V, C>>(
459 op.getmat(), n, w, MILU_VARIANT::ILU);
460 });
461 F::addCreator("ILUn", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
462 const int n = prm.get<int>("ilulevel", 0);
463 const double w = prm.get<double>("relaxation", 1.0);
464 return std::make_shared<ParallelOverlappingILU0<M, V, V, C>>(
465 op.getmat(), n, w, MILU_VARIANT::ILU);
466 });
467 F::addCreator("DILU", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
469 return std::make_shared<MultithreadDILU<M, V, V>>(op.getmat());
470 });
471 F::addCreator("Jac", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
472 const int n = prm.get<int>("repeats", 1);
473 const double w = prm.get<double>("relaxation", 1.0);
474 return getDummyUpdateWrapper<SeqJac<M, V, V>>(op.getmat(), n, w);
475 });
476 F::addCreator("GS", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
477 const int n = prm.get<int>("repeats", 1);
478 const double w = prm.get<double>("relaxation", 1.0);
479 return getDummyUpdateWrapper<SeqGS<M, V, V>>(op.getmat(), n, w);
480 });
481 F::addCreator("SOR", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
482 const int n = prm.get<int>("repeats", 1);
483 const double w = prm.get<double>("relaxation", 1.0);
484 return getDummyUpdateWrapper<SeqSOR<M, V, V>>(op.getmat(), n, w);
485 });
486 F::addCreator("SSOR", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
487 const int n = prm.get<int>("repeats", 1);
488 const double w = prm.get<double>("relaxation", 1.0);
489 return getDummyUpdateWrapper<SeqSSOR<M, V, V>>(op.getmat(), n, w);
490 });
491
492 // Only add AMG preconditioners to the factory if the operator
493 // is an actual matrix operator.
494 if constexpr (std::is_same_v<O, Dune::MatrixAdapter<M, V, V>>) {
495 F::addCreator("amg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
496 const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
497 if (smoother == "ILU0" || smoother == "ParOverILU0") {
498 using Smoother = SeqILU<M, V, V>;
500 } else if (smoother == "Jac") {
501 using Smoother = SeqJac<M, V, V>;
503 } else if (smoother == "GS") {
504 using Smoother = SeqGS<M, V, V>;
506 } else if (smoother == "DILU") {
507 using Smoother = MultithreadDILU<M, V, V>;
509 } else if (smoother == "SOR") {
510 using Smoother = SeqSOR<M, V, V>;
512 } else if (smoother == "SSOR") {
513 using Smoother = SeqSSOR<M, V, V>;
515 } else if (smoother == "ILUn") {
516 using Smoother = SeqILU<M, V, V>;
518 } else {
519 OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + ".");
520 }
521 });
522 F::addCreator("kamg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
523 const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
524 if (smoother == "ILU0" || smoother == "ParOverILU0") {
525 using Smoother = SeqILU<M, V, V>;
527 } else if (smoother == "Jac") {
528 using Smoother = SeqJac<M, V, V>;
530 } else if (smoother == "SOR") {
531 using Smoother = SeqSOR<M, V, V>;
533 } else if (smoother == "GS") {
534 using Smoother = SeqGS<M, V, V>;
536 } else if (smoother == "SSOR") {
537 using Smoother = SeqSSOR<M, V, V>;
539 } else if (smoother == "ILUn") {
540 using Smoother = SeqILU<M, V, V>;
542 } else {
543 OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + ".");
544 }
545 });
546 F::addCreator("famg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
547 if constexpr (std::is_same_v<typename V::field_type, float>) {
548 OPM_THROW(std::logic_error, "famg requires UMFPack which is not available for floats");
549 return nullptr;
550 } else {
552 Dune::Amg::Parameters parms;
553 parms.setNoPreSmoothSteps(1);
554 parms.setNoPostSmoothSteps(1);
556 }
557 });
558
559#if HAVE_AMGX
560 // Only add AMGX for scalar matrices
561 if constexpr (M::block_type::rows == 1 && M::block_type::cols == 1) {
562 F::addCreator("amgx", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
563 auto prm_copy = prm;
564 prm_copy.put("setup_frequency", Opm::Parameters::Get<Opm::Parameters::CprReuseInterval>());
565 return std::make_shared<Amgx::AmgxPreconditioner<M, V, V>>(op.getmat(), prm_copy);
566 });
567 }
568#endif
569
570#if HAVE_HYPRE
571 // Only add Hypre for scalar matrices
572 if constexpr (M::block_type::rows == 1 && M::block_type::cols == 1 &&
573 std::is_same_v<HYPRE_Real, typename V::field_type>) {
574 F::addCreator("hypre", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
575 return std::make_shared<Hypre::HyprePreconditioner<M, V, V>>(op.getmat(), prm);
576 });
577 }
578#endif
579 }
580 if constexpr (std::is_same_v<O, WellModelMatrixAdapter<M, V, V>>) {
581 F::addCreator(
582 "cprw",
583 [](const O& op, const P& prm, const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
584 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
585 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
586 }
587 using Scalar = typename V::field_type;
590 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
591 op, prm, weightsCalculator, pressureIndex);
592 });
593 }
594
595 F::addCreator(
596 "cpr",
597 [](const O& op, const P& prm, const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
598 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
599 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
600 }
601 using Scalar = typename V::field_type;
603 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
604 op, prm, weightsCalculator, pressureIndex);
605 });
606 F::addCreator(
607 "cprt",
608 [](const O& op, const P& prm, const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
609 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
610 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
611 }
612 using Scalar = typename V::field_type;
614 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
615 op, prm, weightsCalculator, pressureIndex);
616 });
617
618#if HAVE_CUDA
619 F::addCreator("GPUILU0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
620 const double w = prm.get<double>("relaxation", 1.0);
621 using field_type = typename V::field_type;
622 using GpuILU0 = typename gpuistl::
623 GpuSeqILU0<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
624 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, GpuILU0>>(
625 std::make_shared<GpuILU0>(op.getmat(), w));
626 });
627
628 F::addCreator("GPUILU0Float", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
629 const double w = prm.get<double>("relaxation", 1.0);
630 using block_type = typename V::block_type;
631 using VTo = Dune::BlockVector<Dune::FieldVector<float, block_type::dimension>>;
632 using matrix_type_to =
633 typename Dune::BCRSMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
634 using GpuILU0 = typename gpuistl::
635 GpuSeqILU0<matrix_type_to, gpuistl::GpuVector<float>, gpuistl::GpuVector<float>>;
638 auto converted = std::make_shared<Converter>(op.getmat());
639 auto adapted = std::make_shared<Adapter>(std::make_shared<GpuILU0>(converted->getConvertedMatrix(), w));
640 converted->setUnderlyingPreconditioner(adapted);
641 return converted;
642 });
643
644 F::addCreator("GPUJAC", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
645 const double w = prm.get<double>("relaxation", 1.0);
646 using field_type = typename V::field_type;
647 using GPUJac =
648 typename gpuistl::GpuJac<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
649 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, GPUJac>>(
650 std::make_shared<GPUJac>(op.getmat(), w));
651 });
652
653 F::addCreator("OPMGPUILU0", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t) {
654 const bool split_matrix = prm.get<bool>("split_matrix", true);
655 const bool tune_gpu_kernels = prm.get<bool>("tune_gpu_kernels", true);
656 const int mixed_precision_scheme = prm.get<int>("mixed_precision_scheme", 0);
657
658 using field_type = typename V::field_type;
659 using GPUILU0 = typename gpuistl::OpmGpuILU0<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
660
661 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, GPUILU0>>(std::make_shared<GPUILU0>(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme));
662 });
663
664 F::addCreator("GPUDILU", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t) {
665 const bool split_matrix = prm.get<bool>("split_matrix", true);
666 const bool tune_gpu_kernels = prm.get<bool>("tune_gpu_kernels", true);
667 const int mixed_precision_scheme = prm.get<int>("mixed_precision_scheme", 0);
668 using field_type = typename V::field_type;
669 using GPUDILU = typename gpuistl::GpuDILU<M, gpuistl::GpuVector<field_type>, gpuistl::GpuVector<field_type>>;
670 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, GPUDILU>>(std::make_shared<GPUDILU>(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme));
671 });
672
673 F::addCreator("GPUDILUFloat", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t) {
674 const bool split_matrix = prm.get<bool>("split_matrix", true);
675 const bool tune_gpu_kernels = prm.get<bool>("tune_gpu_kernels", true);
676 const int mixed_precision_scheme = prm.get<int>("mixed_precision_scheme", 0);
677
678 using block_type = typename V::block_type;
679 using VTo = Dune::BlockVector<Dune::FieldVector<float, block_type::dimension>>;
680 using matrix_type_to = typename Dune::BCRSMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
681 using GpuDILU = typename gpuistl::GpuDILU<matrix_type_to, gpuistl::GpuVector<float>, gpuistl::GpuVector<float>>;
684 auto converted = std::make_shared<Converter>(op.getmat());
685 auto adapted = std::make_shared<Adapter>(std::make_shared<GpuDILU>(converted->getConvertedMatrix(), split_matrix, tune_gpu_kernels, mixed_precision_scheme));
686 converted->setUnderlyingPreconditioner(adapted);
687 return converted;
688 });
689#endif
690 }
691};
692
693template <class Operator, class Comm>
695{
696}
697
698
699template <class Operator, class Comm>
701PreconditionerFactory<Operator, Comm>::instance()
702{
703 static PreconditionerFactory singleton;
704 return singleton;
705}
706
707template <class Operator, class Comm>
709PreconditionerFactory<Operator, Comm>::doCreate(const Operator& op,
710 const PropertyTree& prm,
711 const std::function<Vector()> weightsCalculator,
712 std::size_t pressureIndex)
713{
714 if (!defAdded_) {
715 StandardPreconditioners<Operator, Comm>::add();
716 defAdded_ = true;
717 }
718 const std::string& type = prm.get<std::string>("type", "ParOverILU0");
719 auto it = creators_.find(type);
720 if (it == creators_.end()) {
721 std::ostringstream msg;
722 msg << "Preconditioner type " << type << " is not registered in the factory. Available types are: ";
723 for (const auto& prec : creators_) {
724 msg << prec.first << ' ';
725 }
726 msg << std::endl;
727 OPM_THROW(std::invalid_argument, msg.str());
728 }
729 return it->second(op, prm, weightsCalculator, pressureIndex);
730}
731
732template <class Operator, class Comm>
734PreconditionerFactory<Operator, Comm>::doCreate(const Operator& op,
735 const PropertyTree& prm,
736 const std::function<Vector()> weightsCalculator,
737 std::size_t pressureIndex,
738 const Comm& comm)
739{
740 if (!defAdded_) {
741 StandardPreconditioners<Operator, Comm>::add();
742 defAdded_ = true;
743 }
744 const std::string& type = prm.get<std::string>("type", "ParOverILU0");
745 auto it = parallel_creators_.find(type);
746 if (it == parallel_creators_.end()) {
747 std::ostringstream msg;
748 msg << "Parallel preconditioner type " << type << " is not registered in the factory. Available types are: ";
749 for (const auto& prec : parallel_creators_) {
750 msg << prec.first << ' ';
751 }
752 msg << std::endl;
753 OPM_THROW(std::invalid_argument, msg.str());
754 }
755 return it->second(op, prm, weightsCalculator, pressureIndex, comm);
756}
757
758template <class Operator, class Comm>
759void
760PreconditionerFactory<Operator, Comm>::doAddCreator(const std::string& type, Creator c)
761{
762 creators_[type] = c;
763}
764
765template <class Operator, class Comm>
766void
767PreconditionerFactory<Operator, Comm>::doAddCreator(const std::string& type, ParCreator c)
768{
769 parallel_creators_[type] = c;
770}
771
772template <class Operator, class Comm>
775 const PropertyTree& prm,
776 const std::function<Vector()>& weightsCalculator,
777 std::size_t pressureIndex)
778{
779 return instance().doCreate(op, prm, weightsCalculator, pressureIndex);
780}
781
782template <class Operator, class Comm>
785 const PropertyTree& prm,
786 const std::function<Vector()>& weightsCalculator,
787 const Comm& comm,
788 std::size_t pressureIndex)
789{
790 return instance().doCreate(op, prm, weightsCalculator, pressureIndex, comm);
791}
792
793
794template <class Operator, class Comm>
797 const PropertyTree& prm,
798 const Comm& comm,
799 std::size_t pressureIndex)
800{
801 return instance().doCreate(op, prm, std::function<Vector()>(), pressureIndex, comm);
802}
803
804template <class Operator, class Comm>
805void
807{
808 instance().doAddCreator(type, creator);
809}
810
811template <class Operator, class Comm>
812void
814{
815 instance().doAddCreator(type, creator);
816}
817
818using CommSeq = Dune::Amg::SequentialInformation;
819
820template<class Scalar, int Dim>
821using OpFSeq = Dune::MatrixAdapter<Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, Dim, Dim>>,
822 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
823 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>>;
824template<class Scalar, int Dim>
825using OpBSeq = Dune::MatrixAdapter<Dune::BCRSMatrix<MatrixBlock<Scalar, Dim, Dim>>,
826 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
827 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>>;
828
829template<class Scalar, int Dim, bool overlap>
831 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
832 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>>;
833
834template<class Scalar, int Dim, bool overlap>
836 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
837 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
838 overlap>;
839
840#if HAVE_MPI
841using CommPar = Dune::OwnerOverlapCopyCommunication<int, int>;
842
843template<class Scalar, int Dim>
844using OpFPar = Dune::OverlappingSchwarzOperator<Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, Dim, Dim>>,
845 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
846 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
847 CommPar>;
848
849template<class Scalar, int Dim>
850using OpBPar = Dune::OverlappingSchwarzOperator<Dune::BCRSMatrix<MatrixBlock<Scalar, Dim, Dim>>,
851 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
852 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
853 CommPar>;
854template<class Scalar, int Dim>
856 Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
857 Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
858 CommPar>;
859
860template<class Scalar, int Dim>
862 Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
863 Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
864 CommPar>;
865
866#define INSTANTIATE_PF_PAR(T,Dim) \
867 template class PreconditionerFactory<OpBSeq<T,Dim>, CommPar>; \
868 template class PreconditionerFactory<OpFPar<T,Dim>, CommPar>; \
869 template class PreconditionerFactory<OpBPar<T,Dim>, CommPar>; \
870 template class PreconditionerFactory<OpGLFPar<T,Dim>, CommPar>; \
871 template class PreconditionerFactory<OpGLBPar<T,Dim>, CommPar>; \
872 template class PreconditionerFactory<OpW<T,Dim, false>, CommPar>; \
873 template class PreconditionerFactory<OpWG<T,Dim, true>, CommPar>; \
874 template class PreconditionerFactory<OpBPar<T,Dim>, CommSeq>; \
875 template class PreconditionerFactory<OpGLBPar<T,Dim>, CommSeq>;
876#endif
877
878#define INSTANTIATE_PF_SEQ(T,Dim) \
879 template class PreconditionerFactory<OpFSeq<T,Dim>, CommSeq>; \
880 template class PreconditionerFactory<OpBSeq<T,Dim>, CommSeq>; \
881 template class PreconditionerFactory<OpW<T,Dim, false>, CommSeq>; \
882 template class PreconditionerFactory<OpWG<T,Dim, true>, CommSeq>;
883
884#if HAVE_MPI
885#define INSTANTIATE_PF(T,Dim) \
886 INSTANTIATE_PF_PAR(T,Dim) \
887 INSTANTIATE_PF_SEQ(T,Dim)
888#else
889#define INSTANTIATE_PF(T,Dim) INSTANTIATE_PF_SEQ(T,Dim)
890#endif
891
892} // namespace Opm
The AMG preconditioner.
Parallel algebraic multigrid based on agglomeration.
Definition amgcpr.hh:88
Definition PreconditionerWithUpdate.hpp:43
The OpenMP thread parallelized DILU preconditioner.
Definition DILU.hpp:53
Dune linear operator that assumes ghost rows are ordered after interior rows.
Definition WellOperators.hpp:402
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition ParallelOverlappingILU0.hpp:131
This is an object factory for creating preconditioners.
Definition PreconditionerFactory.hpp:64
static PrecPtr create(const Operator &op, const PropertyTree &prm, const std::function< Vector()> &weightsCalculator={}, std::size_t pressureIndex=std::numeric_limits< std::size_t >::max())
Create a new serial preconditioner and return a pointer to it.
Definition PreconditionerFactory_impl.hpp:774
std::function< PrecPtr(const Operator &, const PropertyTree &, const std::function< Vector()> &, std::size_t)> Creator
The type of creator functions passed to addCreator().
Definition PreconditionerFactory.hpp:75
static void addCreator(const std::string &type, Creator creator)
Add a creator for a serial preconditioner to the PreconditionerFactory.
Definition PreconditionerFactory_impl.hpp:806
std::shared_ptr< Dune::PreconditionerWithUpdate< Vector, Vector > > PrecPtr
The type of pointer returned by create().
Definition PreconditionerFactory.hpp:71
Definition PressureBhpTransferPolicy.hpp:99
Definition PressureTransferPolicy.hpp:55
Hierarchical collection of key/value pairs.
Definition PropertyTree.hpp:39
T get(const std::string &key) const
Retrieve property value given hierarchical property key.
Definition PropertyTree.cpp:59
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:299
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:225
DILU preconditioner on the GPU.
Definition GpuDILU.hpp:47
Jacobi preconditioner on the GPU.
Definition GpuJac.hpp:47
ILU0 preconditioner on the GPU.
Definition OpmGpuILU0.hpp:48
Makes a CUDA preconditioner available to a CPU simulator.
Definition PreconditionerAdapter.hpp:43
Converts the field type (eg.
Definition PreconditionerConvertFieldTypeAdapter.hpp:86
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
MILU_VARIANT
Definition MILU.hpp:34
@ ILU
Do not perform modified ILU.
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Definition PreconditionerFactory.hpp:43
Definition PreconditionerFactory_impl.hpp:67
Definition PreconditionerFactory_impl.hpp:160
static std::size_t interiorIfGhostLast(const Comm &comm)
Helper method to determine if the local partitioning has the K interior cells from [0,...
Definition PreconditionerFactory_impl.hpp:414