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