21#include <opm/simulators/linalg/PreconditionerFactory.hpp>
23#include <opm/common/ErrorMacros.hpp>
25#include <opm/simulators/linalg/matrixblock.hh>
27#include <opm/simulators/linalg/FlexibleSolver.hpp>
28#include <opm/simulators/linalg/ilufirstelement.hh>
29#include <opm/simulators/linalg/OwningBlockPreconditioner.hpp>
30#include <opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp>
31#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
32#include <opm/simulators/linalg/PressureBhpTransferPolicy.hpp>
33#include <opm/simulators/linalg/PressureTransferPolicy.hpp>
34#include <opm/simulators/linalg/PropertyTree.hpp>
35#include <opm/simulators/linalg/WellOperators.hpp>
37#include <dune/istl/owneroverlapcopy.hh>
38#include <dune/istl/preconditioners.hh>
39#include <dune/istl/paamg/amg.hh>
40#include <dune/istl/paamg/kamg.hh>
41#include <dune/istl/paamg/fastamg.hh>
45template<
class Smoother>
50 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
56 smootherArgs.relaxationFactor = prm.get<
double>(
"relaxation", 1.0);
61template<
class M,
class V,
class C>
67 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
70 const int iluwitdh = prm.get<
int>(
"iluwidth", 0);
72 const MILU_VARIANT milu = convertString2Milu(prm.get<std::string>(
"milutype", std::string(
"ilu")));
77 smootherArgs.relaxationFactor = prm.get<
double>(
"relaxation", 1.0);
82template <
class Operator,
class Comm,
class Matrix,
class Vector>
83typename AMGHelper<Operator, Comm, Matrix, Vector>::Criterion
86 Criterion criterion(15, prm.get<
int>(
"coarsenTarget", 1200));
87 criterion.setDefaultValuesIsotropic(2);
88 criterion.setAlpha(prm.get<
double>(
"alpha", 0.33));
89 criterion.setBeta(prm.get<
double>(
"beta", 1
e-5));
90 criterion.setMaxLevel(prm.get<
int>(
"maxlevel", 15));
91 criterion.setSkipIsolated(prm.get<
bool>(
"skip_isolated",
false));
92 criterion.setNoPreSmoothSteps(prm.get<
int>(
"pre_smooth", 1));
93 criterion.setNoPostSmoothSteps(prm.get<
int>(
"post_smooth", 1));
94 criterion.setDebugLevel(prm.get<
int>(
"verbosity", 0));
98 criterion.setAccumulate(
static_cast<Dune::Amg::AccumulationMode
>(prm.get<
int>(
"accumulate", 1)));
99 criterion.setProlongationDampingFactor(prm.get<
double>(
"prolongationdamping", 1.6));
100 criterion.setMaxDistance(prm.get<
int>(
"maxdistance", 2));
101 criterion.setMaxConnectivity(prm.get<
int>(
"maxconnectivity", 15));
102 criterion.setMaxAggregateSize(prm.get<
int>(
"maxaggsize", 6));
103 criterion.setMinAggregateSize(prm.get<
int>(
"minaggsize", 4));
107template <
class Operator,
class Comm,
class Matrix,
class Vector>
108template <
class Smoother>
109typename AMGHelper<Operator, Comm, Matrix, Vector>::PrecPtr
110AMGHelper<Operator,Comm,Matrix,Vector>::
111makeAmgPreconditioner(
const Operator& op,
112 const PropertyTree& prm,
115 auto crit = criterion(prm);
116 auto sargs = AMGSmootherArgsHelper<Smoother>::args(prm);
119 return std::make_shared<Type>(op, crit, sargs,
120 prm.get<
size_t>(
"max_krylov", 1),
121 prm.get<
double>(
"min_reduction", 1e-1));
124 return std::make_shared<Type>(op, crit, sargs);
128template<
class Operator,
class Comm>
133 using namespace Dune;
137 using M =
typename F::Matrix;
138 using V =
typename F::Vector;
140 F::addCreator(
"ILU0", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
141 return createParILU(
op, prm, comm, 0);
143 F::addCreator(
"ParOverILU0", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
144 return createParILU(
op, prm, comm, prm.get<
int>(
"ilulevel", 0));
146 F::addCreator(
"ILUn", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
147 return createParILU(
op, prm, comm, prm.get<
int>(
"ilulevel", 0));
149 F::addCreator(
"Jac", [](
const O&
op,
const P& prm,
const std::function<
V()>&,
150 std::size_t,
const C& comm) {
151 const int n = prm.get<
int>(
"repeats", 1);
152 const double w = prm.get<
double>(
"relaxation", 1.0);
155 F::addCreator(
"GS", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
156 const int n = prm.get<
int>(
"repeats", 1);
157 const double w = prm.get<
double>(
"relaxation", 1.0);
160 F::addCreator(
"SOR", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
161 const int n = prm.get<
int>(
"repeats", 1);
162 const double w = prm.get<
double>(
"relaxation", 1.0);
165 F::addCreator(
"SSOR", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
166 const int n = prm.get<
int>(
"repeats", 1);
167 const double w = prm.get<
double>(
"relaxation", 1.0);
175 if constexpr (std::is_same_v<O, Dune::OverlappingSchwarzOperator<M, V, V, C>>) {
176 F::addCreator(
"amg", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
177 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
178 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
182 return std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(
op,
crit,
sargs, comm);
184 OPM_THROW(std::invalid_argument,
"Properties: No smoother with name " + smoother +
".");
189 F::addCreator(
"cpr", [](
const O&
op,
const P& prm,
const std::function<
V()>
weightsCalculator, std::size_t pressureIndex,
const C& comm) {
191 if (pressureIndex == std::numeric_limits<std::size_t>::max())
193 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
196 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
op, prm,
weightsCalculator, pressureIndex, comm);
198 F::addCreator(
"cprt", [](
const O&
op,
const P& prm,
const std::function<
V()>
weightsCalculator, std::size_t pressureIndex,
const C& comm) {
200 if (pressureIndex == std::numeric_limits<std::size_t>::max())
202 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
205 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
op, prm,
weightsCalculator, pressureIndex, comm);
208 if constexpr (std::is_same_v<O, WellModelGhostLastMatrixAdapter<M, V, V, true>>) {
209 F::addCreator(
"cprw",
210 [](
const O&
op,
const P& prm,
const std::function<
V()>
weightsCalculator, std::size_t pressureIndex,
const C& comm) {
212 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
213 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
216 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
op, prm,
weightsCalculator, pressureIndex, comm);
226 using M =
typename F::Matrix;
227 using V =
typename F::Vector;
229 const double w = prm.get<
double>(
"relaxation", 1.0);
230 const bool redblack = prm.get<
bool>(
"redblack",
false);
235 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, Comm>>(
238 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, Comm>>(
251 const auto&
is = comm.indexSet();
252 for (
const auto&
ind :
is) {
253 if (Comm::OwnerSet::contains(
ind.local().attribute())) {
267template<
class Operator>
272 using namespace Dune;
274 using C = Dune::Amg::SequentialInformation;
276 using M =
typename F::Matrix;
277 using V =
typename F::Vector;
279 F::addCreator(
"ILU0", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
280 const double w = prm.get<
double>(
"relaxation", 1.0);
281 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
284 F::addCreator(
"ParOverILU0", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
285 const double w = prm.get<
double>(
"relaxation", 1.0);
286 const int n = prm.get<
int>(
"ilulevel", 0);
287 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
290 F::addCreator(
"ILUn", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
291 const int n = prm.get<
int>(
"ilulevel", 0);
292 const double w = prm.get<
double>(
"relaxation", 1.0);
293 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
296 F::addCreator(
"Jac", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
297 const int n = prm.get<
int>(
"repeats", 1);
298 const double w = prm.get<
double>(
"relaxation", 1.0);
301 F::addCreator(
"GS", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
302 const int n = prm.get<
int>(
"repeats", 1);
303 const double w = prm.get<
double>(
"relaxation", 1.0);
306 F::addCreator(
"SOR", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
307 const int n = prm.get<
int>(
"repeats", 1);
308 const double w = prm.get<
double>(
"relaxation", 1.0);
311 F::addCreator(
"SSOR", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
312 const int n = prm.get<
int>(
"repeats", 1);
313 const double w = prm.get<
double>(
"relaxation", 1.0);
319 if constexpr (std::is_same_v<O, Dune::MatrixAdapter<M, V, V>>) {
320 F::addCreator(
"amg", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
321 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
322 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
325 }
else if (smoother ==
"Jac") {
328 }
else if (smoother ==
"SOR") {
331 }
else if (smoother ==
"SSOR") {
334 }
else if (smoother ==
"ILUn") {
339 "Properties: No smoother with name " + smoother +
".");
342 F::addCreator(
"kamg", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
343 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
344 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
347 }
else if (smoother ==
"Jac") {
350 }
else if (smoother ==
"SOR") {
356 }
else if (smoother ==
"SSOR") {
359 }
else if (smoother ==
"ILUn") {
364 "Properties: No smoother with name " + smoother +
".");
367 F::addCreator(
"famg", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
369 Dune::Amg::Parameters
parms;
370 parms.setNoPreSmoothSteps(1);
371 parms.setNoPostSmoothSteps(1);
375 if constexpr (std::is_same_v<O, WellModelMatrixAdapter<M, V, V, false>>) {
376 F::addCreator(
"cprw", [](
const O&
op,
const P& prm,
const std::function<
V()>&
weightsCalculator, std::size_t pressureIndex) {
377 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
378 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
381 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
op, prm,
weightsCalculator, pressureIndex);
385 F::addCreator(
"cpr", [](
const O&
op,
const P& prm,
const std::function<
V()>&
weightsCalculator, std::size_t pressureIndex) {
386 if (pressureIndex == std::numeric_limits<std::size_t>::max())
388 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
391 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
op, prm,
weightsCalculator, pressureIndex);
393 F::addCreator(
"cprt", [](
const O&
op,
const P& prm,
const std::function<
V()>&
weightsCalculator, std::size_t pressureIndex) {
394 if (pressureIndex == std::numeric_limits<std::size_t>::max())
396 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
399 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
op, prm,
weightsCalculator, pressureIndex);
404template <
class Operator,
class Comm>
410template <
class Operator,
class Comm>
411PreconditionerFactory<Operator,Comm>&
412PreconditionerFactory<Operator,Comm>::instance()
414 static PreconditionerFactory singleton;
418template <
class Operator,
class Comm>
420PreconditionerFactory<Operator,Comm>::
421doCreate(
const Operator& op,
const PropertyTree& prm,
422 const std::function<Vector()> weightsCalculator,
423 std::size_t pressureIndex)
426 StandardPreconditioners<Operator,Comm>::add();
429 const std::string& type = prm.get<std::string>(
"type",
"ParOverILU0");
430 auto it = creators_.find(type);
431 if (it == creators_.end()) {
432 std::ostringstream msg;
433 msg <<
"Preconditioner type " << type <<
" is not registered in the factory. Available types are: ";
434 for (
const auto& prec : creators_) {
435 msg << prec.first <<
' ';
438 OPM_THROW(std::invalid_argument, msg.str());
440 return it->second(op, prm, weightsCalculator, pressureIndex);
443template <
class Operator,
class Comm>
445PreconditionerFactory<Operator,Comm>::
446doCreate(
const Operator& op,
const PropertyTree& prm,
447 const std::function<Vector()> weightsCalculator,
448 std::size_t pressureIndex,
const Comm& comm)
451 StandardPreconditioners<Operator,Comm>::add();
454 const std::string& type = prm.get<std::string>(
"type",
"ParOverILU0");
455 auto it = parallel_creators_.find(type);
456 if (it == parallel_creators_.end()) {
457 std::ostringstream msg;
458 msg <<
"Parallel preconditioner type " << type
459 <<
" is not registered in the factory. Available types are: ";
460 for (
const auto& prec : parallel_creators_) {
461 msg << prec.first <<
' ';
464 OPM_THROW(std::invalid_argument, msg.str());
466 return it->second(op, prm, weightsCalculator, pressureIndex, comm);
469template <
class Operator,
class Comm>
470void PreconditionerFactory<Operator,Comm>::
471doAddCreator(
const std::string& type, Creator c)
476template <
class Operator,
class Comm>
477void PreconditionerFactory<Operator,Comm>::
478doAddCreator(
const std::string& type, ParCreator c)
480 parallel_creators_[type] = c;
483template <
class Operator,
class Comm>
488 std::size_t pressureIndex)
493template <
class Operator,
class Comm>
498 std::size_t pressureIndex)
504template <
class Operator,
class Comm>
508 std::size_t pressureIndex)
510 return instance().doCreate(
op, prm, std::function<Vector()>(), pressureIndex, comm);
513template <
class Operator,
class Comm>
517 instance().doAddCreator(type,
creator);
520template <
class Operator,
class Comm>
524 instance().doAddCreator(type,
creator);
527using CommSeq = Dune::Amg::SequentialInformation;
530using OpFSeq = Dune::MatrixAdapter<Dune::BCRSMatrix<Dune::FieldMatrix<double,Dim,Dim>>,
531 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
532 Dune::BlockVector<Dune::FieldVector<double,Dim>>>;
534using OpBSeq = Dune::MatrixAdapter<Dune::BCRSMatrix<Opm::MatrixBlock<double,Dim,Dim>>,
535 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
536 Dune::BlockVector<Dune::FieldVector<double,Dim>>>;
538template<
int Dim,
bool overlap>
540 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
541 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
544template<
int Dim,
bool overlap>
546 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
547 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
551using CommPar = Dune::OwnerOverlapCopyCommunication<int,int>;
554using OpFPar = Dune::OverlappingSchwarzOperator<Dune::BCRSMatrix<Dune::FieldMatrix<double,Dim,Dim>>,
555 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
556 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
560using OpBPar = Dune::OverlappingSchwarzOperator<Dune::BCRSMatrix<MatrixBlock<double,Dim,Dim>>,
561 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
562 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
565#define INSTANCE_PF_PAR(Dim) \
566 template class PreconditionerFactory<OpBSeq<Dim>,CommPar>; \
567 template class PreconditionerFactory<OpFPar<Dim>,CommPar>; \
568 template class PreconditionerFactory<OpBPar<Dim>,CommPar>; \
569 template class PreconditionerFactory<OpW<Dim,false>,CommPar>; \
570 template class PreconditionerFactory<OpWG<Dim,true>,CommPar>; \
571 template class PreconditionerFactory<OpBPar<Dim>,CommSeq>;
574#define INSTANCE_PF_SEQ(Dim) \
575 template class PreconditionerFactory<OpFSeq<Dim>,CommSeq>; \
576 template class PreconditionerFactory<OpBSeq<Dim>,CommSeq>; \
577 template class PreconditionerFactory<OpW<Dim,false>,CommSeq>; \
578 template class PreconditionerFactory<OpWG<Dim,true>,CommSeq>;
581#define INSTANCE_PF(Dim) \
582 INSTANCE_PF_PAR(Dim) \
585#define INSTANCE_PF(Dim) \
Parallel algebraic multigrid based on agglomeration.
Definition amgcpr.hh:87
Definition PreconditionerWithUpdate.hpp:40
Definition AquiferInterface.hpp:35
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:486
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:515
std::shared_ptr< Dune::PreconditionerWithUpdate< Vector, Vector > > PrecPtr
The type of pointer returned by create().
Definition PreconditionerFactory.hpp:71
Definition PropertyTree.hpp:37
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
MILU_VARIANT
Definition MILU.hpp:34
@ ILU
Do not perform modified ILU.
Definition PreconditionerFactory.hpp:43
Definition PreconditionerFactory_impl.hpp:47
Definition PreconditionerFactory_impl.hpp:130
static 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:247