22 #ifndef OPM_PRECONDITIONERFACTORY_HEADER
23 #define OPM_PRECONDITIONERFACTORY_HEADER
25 #include <opm/simulators/linalg/OwningBlockPreconditioner.hpp>
26 #include <opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp>
27 #include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
28 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
29 #include <opm/simulators/linalg/PropertyTree.hpp>
32 #include <dune/istl/paamg/amg.hh>
33 #include <dune/istl/paamg/kamg.hh>
34 #include <dune/istl/paamg/fastamg.hh>
35 #include <dune/istl/preconditioners.hh>
49 template <
class Operator,
class Comm>
54 using Matrix =
typename Operator::matrix_type;
55 using Vector =
typename Operator::domain_type;
58 using PrecPtr = std::shared_ptr<Dune::PreconditionerWithUpdate<Vector, Vector>>;
62 const std::function<Vector()>&, std::size_t)>;
64 const std::function<Vector()>&, std::size_t,
const Comm&)>;
72 const std::function<Vector()>& weightsCalculator = {},
73 std::size_t pressureIndex = std::numeric_limits<std::size_t>::max())
75 return instance().doCreate(op, prm, weightsCalculator, pressureIndex);
85 const std::function<Vector()>& weightsCalculator,
const Comm& comm,
86 std::size_t pressureIndex = std::numeric_limits<std::size_t>::max())
88 return instance().doCreate(op, prm, weightsCalculator, pressureIndex, comm);
97 std::size_t pressureIndex = std::numeric_limits<std::size_t>::max())
99 return instance().doCreate(op, prm, std::function<Vector()>(), pressureIndex, comm);
110 instance().doAddCreator(type, creator);
120 static void addCreator(
const std::string& type, ParCreator creator)
122 instance().doAddCreator(type, creator);
126 = Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricDependency<Matrix, Dune::Amg::FirstDiagonal>>;
127 using Criterion = Dune::Amg::CoarsenCriterion<CriterionBase>;
132 Criterion criterion(15, prm.get<
int>(
"coarsenTarget", 1200));
133 criterion.setDefaultValuesIsotropic(2);
134 criterion.setAlpha(prm.get<
double>(
"alpha", 0.33));
135 criterion.setBeta(prm.get<
double>(
"beta", 1e-5));
136 criterion.setMaxLevel(prm.get<
int>(
"maxlevel", 15));
137 criterion.setSkipIsolated(prm.get<
bool>(
"skip_isolated",
false));
138 criterion.setNoPreSmoothSteps(prm.get<
int>(
"pre_smooth", 1));
139 criterion.setNoPostSmoothSteps(prm.get<
int>(
"post_smooth", 1));
140 criterion.setDebugLevel(prm.get<
int>(
"verbosity", 0));
144 criterion.setAccumulate(
static_cast<Dune::Amg::AccumulationMode
>(prm.get<
int>(
"accumulate", 1)));
145 criterion.setProlongationDampingFactor(prm.get<
double>(
"prolongationdamping", 1.6));
146 criterion.setMaxDistance(prm.get<
int>(
"maxdistance", 2));
147 criterion.setMaxConnectivity(prm.get<
int>(
"maxconnectivity", 15));
148 criterion.setMaxAggregateSize(prm.get<
int>(
"maxaggsize", 6));
149 criterion.setMinAggregateSize(prm.get<
int>(
"minaggsize", 4));
156 template <
typename X>
struct Id {
using Type = X; };
158 template <
typename Smoother>
159 static auto amgSmootherArgs(
const PropertyTree& prm)
161 return amgSmootherArgs(prm, Id<Smoother>());
164 template <
typename Smoother>
165 static auto amgSmootherArgs(
const PropertyTree& prm,
168 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
169 SmootherArgs smootherArgs;
170 smootherArgs.iterations = prm.get<
int>(
"iterations", 1);
174 smootherArgs.relaxationFactor = prm.get<
double>(
"relaxation", 1.0);
178 static auto amgSmootherArgs(
const PropertyTree& prm,
182 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
183 SmootherArgs smootherArgs;
184 smootherArgs.iterations = prm.get<
int>(
"iterations", 1);
185 const int iluwitdh = prm.get<
int>(
"iluwidth", 0);
186 smootherArgs.setN(iluwitdh);
187 const MILU_VARIANT milu = convertString2Milu(prm.get<std::string>(
"milutype", std::string(
"ilu")));
188 smootherArgs.setMilu(milu);
192 smootherArgs.relaxationFactor = prm.get<
double>(
"relaxation", 1.0);
196 template <
class Smoother>
197 static PrecPtr makeAmgPreconditioner(
const Operator& op,
const PropertyTree& prm,
bool useKamg =
false)
199 auto crit = amgCriterion(prm);
200 auto sargs = amgSmootherArgs<Smoother>(prm);
202 return std::make_shared<
207 prm.get<
size_t>(
"max_krylov", 1),
208 prm.get<
double>(
"min_reduction", 1e-1) );
210 return std::make_shared<Dune::Amg::AMGCPR<Operator, Vector, Smoother>>(op, crit, sargs);
218 template <
class CommArg>
219 static size_t interiorIfGhostLast(
const CommArg& comm)
221 size_t interior_count = 0;
222 size_t highest_interior_index = 0;
223 const auto& is = comm.indexSet();
224 for (
const auto& ind : is) {
225 if (CommArg::OwnerSet::contains(ind.local().attribute())) {
227 highest_interior_index = std::max(highest_interior_index, ind.local().local());
230 if (highest_interior_index + 1 == interior_count) {
231 return interior_count;
238 createParILU(
const Operator& op,
const PropertyTree& prm,
const Comm& comm,
const int ilulevel)
240 const double w = prm.get<
double>(
"relaxation", 1.0);
241 const bool redblack = prm.get<
bool>(
"redblack",
false);
242 const bool reorder_spheres = prm.get<
bool>(
"reorder_spheres",
false);
245 const size_t num_interior = interiorIfGhostLast(comm);
246 return std::make_shared<Opm::ParallelOverlappingILU0<Matrix, Vector, Vector, Comm>>(
249 return std::make_shared<Opm::ParallelOverlappingILU0<Matrix, Vector, Vector, Comm>>(
257 template <
class CommArg>
258 void addStandardPreconditioners(
const CommArg*)
260 using namespace Dune;
264 using P = PropertyTree;
266 doAddCreator(
"ILU0", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t,
const C& comm) {
267 return createParILU(op, prm, comm, 0);
269 doAddCreator(
"ParOverILU0", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t,
const C& comm) {
270 return createParILU(op, prm, comm, prm.get<
int>(
"ilulevel", 0));
272 doAddCreator(
"ILUn", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t,
const C& comm) {
273 return createParILU(op, prm, comm, prm.get<
int>(
"ilulevel", 0));
275 doAddCreator(
"Jac", [](
const O& op,
const P& prm,
const std::function<Vector()>&,
276 std::size_t,
const C& comm) {
277 const int n = prm.get<
int>(
"repeats", 1);
278 const double w = prm.get<
double>(
"relaxation", 1.0);
279 return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqJac<M, V, V>>>(comm, op.getmat(), n, w);
281 doAddCreator(
"GS", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t,
const C& comm) {
282 const int n = prm.get<
int>(
"repeats", 1);
283 const double w = prm.get<
double>(
"relaxation", 1.0);
284 return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqGS<M, V, V>>>(comm, op.getmat(), n, w);
286 doAddCreator(
"SOR", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t,
const C& comm) {
287 const int n = prm.get<
int>(
"repeats", 1);
288 const double w = prm.get<
double>(
"relaxation", 1.0);
289 return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSOR<M, V, V>>>(comm, op.getmat(), n, w);
291 doAddCreator(
"SSOR", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t,
const C& comm) {
292 const int n = prm.get<
int>(
"repeats", 1);
293 const double w = prm.get<
double>(
"relaxation", 1.0);
294 return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSSOR<M, V, V>>>(comm, op.getmat(), n, w);
301 if constexpr (std::is_same_v<O, Dune::OverlappingSchwarzOperator<M, V, V, C>>) {
302 doAddCreator(
"amg", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t,
const C& comm) {
303 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
304 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
306 auto crit = amgCriterion(prm);
307 auto sargs = amgSmootherArgs<Smoother>(prm);
308 return std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
310 OPM_THROW(std::invalid_argument,
"Properties: No smoother with name " << smoother <<
".");
315 doAddCreator(
"cpr", [](
const O& op,
const P& prm,
const std::function<Vector()> weightsCalculator, std::size_t pressureIndex,
const C& comm) {
316 assert(weightsCalculator);
317 if (pressureIndex == std::numeric_limits<std::size_t>::max())
319 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
321 return std::make_shared<OwningTwoLevelPreconditioner<O, V, false, Comm>>(op, prm, weightsCalculator, pressureIndex, comm);
323 doAddCreator(
"cprt", [](
const O& op,
const P& prm,
const std::function<Vector()> weightsCalculator, std::size_t pressureIndex,
const C& comm) {
324 assert(weightsCalculator);
325 if (pressureIndex == std::numeric_limits<std::size_t>::max())
327 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
329 return std::make_shared<OwningTwoLevelPreconditioner<O, V, true, Comm>>(op, prm, weightsCalculator, pressureIndex, comm);
335 void addStandardPreconditioners(
const Dune::Amg::SequentialInformation*)
337 using namespace Dune;
341 using P = PropertyTree;
342 doAddCreator(
"ILU0", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t) {
343 const double w = prm.get<
double>(
"relaxation", 1.0);
344 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
347 doAddCreator(
"ParOverILU0", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t) {
348 const double w = prm.get<
double>(
"relaxation", 1.0);
349 const int n = prm.get<
int>(
"ilulevel", 0);
350 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
353 doAddCreator(
"ILUn", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t) {
354 const int n = prm.get<
int>(
"ilulevel", 0);
355 const double w = prm.get<
double>(
"relaxation", 1.0);
356 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
359 doAddCreator(
"Jac", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t) {
360 const int n = prm.get<
int>(
"repeats", 1);
361 const double w = prm.get<
double>(
"relaxation", 1.0);
362 return wrapPreconditioner<SeqJac<M, V, V>>(op.getmat(), n, w);
364 doAddCreator(
"GS", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t) {
365 const int n = prm.get<
int>(
"repeats", 1);
366 const double w = prm.get<
double>(
"relaxation", 1.0);
367 return wrapPreconditioner<SeqGS<M, V, V>>(op.getmat(), n, w);
369 doAddCreator(
"SOR", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t) {
370 const int n = prm.get<
int>(
"repeats", 1);
371 const double w = prm.get<
double>(
"relaxation", 1.0);
372 return wrapPreconditioner<SeqSOR<M, V, V>>(op.getmat(), n, w);
374 doAddCreator(
"SSOR", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t) {
375 const int n = prm.get<
int>(
"repeats", 1);
376 const double w = prm.get<
double>(
"relaxation", 1.0);
377 return wrapPreconditioner<SeqSSOR<M, V, V>>(op.getmat(), n, w);
382 if constexpr (std::is_same_v<O, Dune::MatrixAdapter<M, V, V>>) {
383 doAddCreator(
"amg", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t) {
384 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
385 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
386 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
387 using Smoother = SeqILU<M, V, V>;
389 using Smoother = SeqILU0<M, V, V>;
391 return makeAmgPreconditioner<Smoother>(op, prm);
392 }
else if (smoother ==
"Jac") {
393 using Smoother = SeqJac<M, V, V>;
394 return makeAmgPreconditioner<Smoother>(op, prm);
395 }
else if (smoother ==
"SOR") {
396 using Smoother = SeqSOR<M, V, V>;
397 return makeAmgPreconditioner<Smoother>(op, prm);
398 }
else if (smoother ==
"SSOR") {
399 using Smoother = SeqSSOR<M, V, V>;
400 return makeAmgPreconditioner<Smoother>(op, prm);
401 }
else if (smoother ==
"ILUn") {
402 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
403 using Smoother = SeqILU<M, V, V>;
405 using Smoother = SeqILUn<M, V, V>;
407 return makeAmgPreconditioner<Smoother>(op, prm);
409 OPM_THROW(std::invalid_argument,
"Properties: No smoother with name " << smoother <<
".");
412 doAddCreator(
"kamg", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t) {
413 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
414 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
415 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
416 using Smoother = SeqILU<M, V, V>;
418 using Smoother = SeqILU0<M, V, V>;
420 return makeAmgPreconditioner<Smoother>(op, prm,
true);
421 }
else if (smoother ==
"Jac") {
422 using Smoother = SeqJac<M, V, V>;
423 return makeAmgPreconditioner<Smoother>(op, prm,
true);
424 }
else if (smoother ==
"SOR") {
425 using Smoother = SeqSOR<M, V, V>;
426 return makeAmgPreconditioner<Smoother>(op, prm,
true);
430 }
else if (smoother ==
"SSOR") {
431 using Smoother = SeqSSOR<M, V, V>;
432 return makeAmgPreconditioner<Smoother>(op, prm,
true);
433 }
else if (smoother ==
"ILUn") {
434 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
435 using Smoother = SeqILU<M, V, V>;
437 using Smoother = SeqILUn<M, V, V>;
439 return makeAmgPreconditioner<Smoother>(op, prm,
true);
441 OPM_THROW(std::invalid_argument,
"Properties: No smoother with name " << smoother <<
".");
444 doAddCreator(
"famg", [](
const O& op,
const P& prm,
const std::function<Vector()>&, std::size_t) {
445 auto crit = amgCriterion(prm);
446 Dune::Amg::Parameters parms;
447 parms.setNoPreSmoothSteps(1);
448 parms.setNoPostSmoothSteps(1);
449 return wrapPreconditioner<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
452 doAddCreator(
"cpr", [](
const O& op,
const P& prm,
const std::function<Vector()>& weightsCalculator, std::size_t pressureIndex) {
453 if (pressureIndex == std::numeric_limits<std::size_t>::max())
455 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
457 return std::make_shared<OwningTwoLevelPreconditioner<O, V, false>>(op, prm, weightsCalculator, pressureIndex);
459 doAddCreator(
"cprt", [](
const O& op,
const P& prm,
const std::function<Vector()>& weightsCalculator, std::size_t pressureIndex) {
460 if (pressureIndex == std::numeric_limits<std::size_t>::max())
462 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
464 return std::make_shared<OwningTwoLevelPreconditioner<O, V, true>>(op, prm, weightsCalculator, pressureIndex);
471 static PreconditionerFactory& instance()
473 static PreconditionerFactory singleton;
478 PreconditionerFactory()
480 Comm* dummy =
nullptr;
481 addStandardPreconditioners(dummy);
485 PrecPtr doCreate(
const Operator& op,
const PropertyTree& prm,
486 const std::function<Vector()> weightsCalculator,
487 std::size_t pressureIndex)
489 const std::string& type = prm.get<std::string>(
"type",
"ParOverILU0");
490 auto it = creators_.find(type);
491 if (it == creators_.end()) {
492 std::ostringstream msg;
493 msg <<
"Preconditioner type " << type <<
" is not registered in the factory. Available types are: ";
494 for (
const auto& prec : creators_) {
495 msg << prec.first <<
' ';
498 OPM_THROW(std::invalid_argument, msg.str());
500 return it->second(op, prm, weightsCalculator, pressureIndex);
503 PrecPtr doCreate(
const Operator& op,
const PropertyTree& prm,
504 const std::function<Vector()> weightsCalculator,
505 std::size_t pressureIndex,
const Comm& comm)
507 const std::string& type = prm.get<std::string>(
"type",
"ParOverILU0");
508 auto it = parallel_creators_.find(type);
509 if (it == parallel_creators_.end()) {
510 std::ostringstream msg;
511 msg <<
"Parallel preconditioner type " << type
512 <<
" is not registered in the factory. Available types are: ";
513 for (
const auto& prec : parallel_creators_) {
514 msg << prec.first <<
' ';
517 OPM_THROW(std::invalid_argument, msg.str());
519 return it->second(op, prm, weightsCalculator, pressureIndex, comm);
523 void doAddCreator(
const std::string& type,
Creator c)
529 void doAddCreator(
const std::string& type, ParCreator c)
531 parallel_creators_[type] = c;
535 std::map<std::string, Creator> creators_;
536 std::map<std::string, ParCreator> parallel_creators_;
Definition: PreconditionerWithUpdate.hpp:40
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition: ParallelOverlappingILU0.hpp:611
This is an object factory for creating preconditioners.
Definition: PreconditionerFactory.hpp:51
static PrecPtr create(const Operator &op, const PropertyTree &prm, const std::function< Vector()> &weightsCalculator, const Comm &comm, std::size_t pressureIndex=std::numeric_limits< std::size_t >::max())
Create a new parallel preconditioner and return a pointer to it.
Definition: PreconditionerFactory.hpp:84
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:62
static void addCreator(const std::string &type, ParCreator creator)
Add a creator for a parallel preconditioner to the PreconditionerFactory.
Definition: PreconditionerFactory.hpp:120
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.hpp:71
std::shared_ptr< Dune::PreconditionerWithUpdate< Vector, Vector > > PrecPtr
The type of pointer returned by create().
Definition: PreconditionerFactory.hpp:58
static PrecPtr create(const Operator &op, const PropertyTree &prm, const Comm &comm, std::size_t pressureIndex=std::numeric_limits< std::size_t >::max())
Create a new parallel preconditioner and return a pointer to it.
Definition: PreconditionerFactory.hpp:96
static void addCreator(const std::string &type, Creator creator)
Add a creator for a serial preconditioner to the PreconditionerFactory.
Definition: PreconditionerFactory.hpp:108
typename Operator::matrix_type Matrix
Linear algebra types.
Definition: PreconditionerFactory.hpp:54
Definition: PropertyTree.hpp:37
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
MILU_VARIANT
Definition: ParallelOverlappingILU0.hpp:47
@ ILU
Do not perform modified ILU.