20#include <opm/common/ErrorMacros.hpp>
21#include <opm/common/TimingMacros.hpp>
23#include <opm/simulators/linalg/PreconditionerFactory.hpp>
26#include <opm/simulators/linalg/matrixblock.hh>
28#include <opm/simulators/linalg/FlexibleSolver.hpp>
29#include <opm/simulators/linalg/ilufirstelement.hh>
30#include <opm/simulators/linalg/OwningBlockPreconditioner.hpp>
31#include <opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp>
32#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
33#include <opm/simulators/linalg/PressureBhpTransferPolicy.hpp>
34#include <opm/simulators/linalg/PressureTransferPolicy.hpp>
35#include <opm/simulators/linalg/PropertyTree.hpp>
36#include <opm/simulators/linalg/WellOperators.hpp>
38#include <dune/istl/owneroverlapcopy.hh>
39#include <dune/istl/preconditioners.hh>
40#include <dune/istl/paamg/amg.hh>
41#include <dune/istl/paamg/kamg.hh>
42#include <dune/istl/paamg/fastamg.hh>
46#include <opm/simulators/linalg/cuistl/CuSeqILU0.hpp>
47#include <opm/simulators/linalg/cuistl/PreconditionerAdapter.hpp>
48#include <opm/simulators/linalg/cuistl/PreconditionerConvertFieldTypeAdapter.hpp>
49#include <opm/simulators/linalg/cuistl/CuBlockPreconditioner.hpp>
56template<
class Smoother>
61 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
67 smootherArgs.relaxationFactor = prm.get<
double>(
"relaxation", 1.0);
72template<
class M,
class V,
class C>
78 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
81 const int iluwitdh = prm.get<
int>(
"iluwidth", 0);
83 const MILU_VARIANT milu = convertString2Milu(prm.get<std::string>(
"milutype", std::string(
"ilu")));
88 smootherArgs.relaxationFactor = prm.get<
double>(
"relaxation", 1.0);
93template <
class Operator,
class Comm,
class Matrix,
class Vector>
94typename AMGHelper<Operator, Comm, Matrix, Vector>::Criterion
97 Criterion criterion(15, prm.get<
int>(
"coarsenTarget", 1200));
98 criterion.setDefaultValuesIsotropic(2);
99 criterion.setAlpha(prm.get<
double>(
"alpha", 0.33));
100 criterion.setBeta(prm.get<
double>(
"beta", 1
e-5));
101 criterion.setMaxLevel(prm.get<
int>(
"maxlevel", 15));
102 criterion.setSkipIsolated(prm.get<
bool>(
"skip_isolated",
false));
103 criterion.setNoPreSmoothSteps(prm.get<
int>(
"pre_smooth", 1));
104 criterion.setNoPostSmoothSteps(prm.get<
int>(
"post_smooth", 1));
105 criterion.setDebugLevel(prm.get<
int>(
"verbosity", 0));
109 criterion.setAccumulate(
static_cast<Dune::Amg::AccumulationMode
>(prm.get<
int>(
"accumulate", 1)));
110 criterion.setProlongationDampingFactor(prm.get<
double>(
"prolongationdamping", 1.6));
111 criterion.setMaxDistance(prm.get<
int>(
"maxdistance", 2));
112 criterion.setMaxConnectivity(prm.get<
int>(
"maxconnectivity", 15));
113 criterion.setMaxAggregateSize(prm.get<
int>(
"maxaggsize", 6));
114 criterion.setMinAggregateSize(prm.get<
int>(
"minaggsize", 4));
118template <
class Operator,
class Comm,
class Matrix,
class Vector>
119template <
class Smoother>
120typename AMGHelper<Operator, Comm, Matrix, Vector>::PrecPtr
121AMGHelper<Operator,Comm,Matrix,Vector>::
122makeAmgPreconditioner(
const Operator& op,
123 const PropertyTree& prm,
126 auto crit = criterion(prm);
127 auto sargs = AMGSmootherArgsHelper<Smoother>::args(prm);
130 return std::make_shared<Type>(op, crit, sargs,
131 prm.get<std::size_t>(
"max_krylov", 1),
132 prm.get<
double>(
"min_reduction", 1e-1));
135 return std::make_shared<Type>(op, crit, sargs);
139template<
class Operator,
class Comm>
144 using namespace Dune;
148 using M =
typename F::Matrix;
149 using V =
typename F::Vector;
151 F::addCreator(
"ILU0", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
152 return createParILU(
op, prm, comm, 0);
154 F::addCreator(
"ParOverILU0", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
155 return createParILU(
op, prm, comm, prm.get<
int>(
"ilulevel", 0));
157 F::addCreator(
"ILUn", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
158 return createParILU(
op, prm, comm, prm.get<
int>(
"ilulevel", 0));
160 F::addCreator(
"Jac", [](
const O&
op,
const P& prm,
const std::function<
V()>&,
161 std::size_t,
const C& comm) {
162 const int n = prm.get<
int>(
"repeats", 1);
163 const double w = prm.get<
double>(
"relaxation", 1.0);
166 F::addCreator(
"GS", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
167 const int n = prm.get<
int>(
"repeats", 1);
168 const double w = prm.get<
double>(
"relaxation", 1.0);
171 F::addCreator(
"SOR", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
172 const int n = prm.get<
int>(
"repeats", 1);
173 const double w = prm.get<
double>(
"relaxation", 1.0);
176 F::addCreator(
"SSOR", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
177 const int n = prm.get<
int>(
"repeats", 1);
178 const double w = prm.get<
double>(
"relaxation", 1.0);
186 if constexpr (std::is_same_v<O, Dune::OverlappingSchwarzOperator<M, V, V, C>>) {
187 F::addCreator(
"amg", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
188 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
189 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
193 return std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(
op,
crit,
sargs, comm);
195 OPM_THROW(std::invalid_argument,
"Properties: No smoother with name " + smoother +
".");
200 F::addCreator(
"cpr", [](
const O&
op,
const P& prm,
const std::function<
V()>
weightsCalculator, std::size_t pressureIndex,
const C& comm) {
202 if (pressureIndex == std::numeric_limits<std::size_t>::max())
204 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
207 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
op, prm,
weightsCalculator, pressureIndex, comm);
209 F::addCreator(
"cprt", [](
const O&
op,
const P& prm,
const std::function<
V()>
weightsCalculator, std::size_t pressureIndex,
const C& comm) {
211 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);
219 if constexpr (std::is_same_v<O, WellModelGhostLastMatrixAdapter<M, V, V, true>>) {
220 F::addCreator(
"cprw",
221 [](
const O&
op,
const P& prm,
const std::function<
V()>
weightsCalculator, std::size_t pressureIndex,
const C& comm) {
223 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
224 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
227 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
op, prm,
weightsCalculator, pressureIndex, comm);
232 F::addCreator(
"CUILU0", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t,
const C& comm) {
233 const double w = prm.get<
double>(
"relaxation", 1.0);
234 using field_type =
typename V::field_type;
236 auto cuILU0 = std::make_shared<CuILU0>(
op.getmat(),
w);
238 auto adapted = std::make_shared<Opm::cuistl::PreconditionerAdapter<V, V, CuILU0>>(
cuILU0);
239 auto wrapped = std::make_shared<Opm::cuistl::CuBlockPreconditioner<V, V, Comm>>(
adapted, comm);
250 using M =
typename F::Matrix;
251 using V =
typename F::Vector;
253 const double w = prm.get<
double>(
"relaxation", 1.0);
254 const bool redblack = prm.get<
bool>(
"redblack",
false);
259 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, Comm>>(
262 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, Comm>>(
275 const auto&
is = comm.indexSet();
276 for (
const auto&
ind :
is) {
277 if (Comm::OwnerSet::contains(
ind.local().attribute())) {
291template<
class Operator>
296 using namespace Dune;
298 using C = Dune::Amg::SequentialInformation;
300 using M =
typename F::Matrix;
301 using V =
typename F::Vector;
303 F::addCreator(
"ILU0", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
304 const double w = prm.get<
double>(
"relaxation", 1.0);
305 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
308 F::addCreator(
"ParOverILU0", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
309 const double w = prm.get<
double>(
"relaxation", 1.0);
310 const int n = prm.get<
int>(
"ilulevel", 0);
311 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
314 F::addCreator(
"ILUn", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
315 const int n = prm.get<
int>(
"ilulevel", 0);
316 const double w = prm.get<
double>(
"relaxation", 1.0);
317 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
320 F::addCreator(
"Jac", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
321 const int n = prm.get<
int>(
"repeats", 1);
322 const double w = prm.get<
double>(
"relaxation", 1.0);
325 F::addCreator(
"GS", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
326 const int n = prm.get<
int>(
"repeats", 1);
327 const double w = prm.get<
double>(
"relaxation", 1.0);
330 F::addCreator(
"SOR", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
331 const int n = prm.get<
int>(
"repeats", 1);
332 const double w = prm.get<
double>(
"relaxation", 1.0);
335 F::addCreator(
"SSOR", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
336 const int n = prm.get<
int>(
"repeats", 1);
337 const double w = prm.get<
double>(
"relaxation", 1.0);
343 if constexpr (std::is_same_v<O, Dune::MatrixAdapter<M, V, V>>) {
344 F::addCreator(
"amg", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
345 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
346 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
349 }
else if (smoother ==
"Jac") {
352 }
else if (smoother ==
"SOR") {
355 }
else if (smoother ==
"SSOR") {
358 }
else if (smoother ==
"ILUn") {
363 "Properties: No smoother with name " + smoother +
".");
366 F::addCreator(
"kamg", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
367 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
368 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
371 }
else if (smoother ==
"Jac") {
374 }
else if (smoother ==
"SOR") {
380 }
else if (smoother ==
"SSOR") {
383 }
else if (smoother ==
"ILUn") {
388 "Properties: No smoother with name " + smoother +
".");
391 F::addCreator(
"famg", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
393 Dune::Amg::Parameters
parms;
394 parms.setNoPreSmoothSteps(1);
395 parms.setNoPostSmoothSteps(1);
399 if constexpr (std::is_same_v<O, WellModelMatrixAdapter<M, V, V, false>>) {
400 F::addCreator(
"cprw", [](
const O&
op,
const P& prm,
const std::function<
V()>&
weightsCalculator, std::size_t pressureIndex) {
401 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
402 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
405 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
op, prm,
weightsCalculator, pressureIndex);
409 F::addCreator(
"cpr", [](
const O&
op,
const P& prm,
const std::function<
V()>&
weightsCalculator, std::size_t pressureIndex) {
410 if (pressureIndex == std::numeric_limits<std::size_t>::max())
412 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
415 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
op, prm,
weightsCalculator, pressureIndex);
417 F::addCreator(
"cprt", [](
const O&
op,
const P& prm,
const std::function<
V()>&
weightsCalculator, std::size_t pressureIndex) {
418 if (pressureIndex == std::numeric_limits<std::size_t>::max())
420 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
423 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
op, prm,
weightsCalculator, pressureIndex);
427 F::addCreator(
"CUILU0", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
428 const double w = prm.get<
double>(
"relaxation", 1.0);
429 using field_type =
typename V::field_type;
431 return std::make_shared<Opm::cuistl::PreconditionerAdapter<V, V, CuILU0>>(std::make_shared<CuILU0>(
op.getmat(),
w));
434 F::addCreator(
"CUILU0Float", [](
const O&
op,
const P& prm,
const std::function<
V()>&, std::size_t) {
435 const double w = prm.get<
double>(
"relaxation", 1.0);
436 using block_type =
typename V::block_type;
437 using VTo = Dune::BlockVector<Dune::FieldVector<float, block_type::dimension>>;
438 using matrix_type_to =
typename Dune::BCRSMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
442 auto converted = std::make_shared<Converter>(
op.getmat());
443 auto adapted = std::make_shared<Adapter>(std::make_shared<CuILU0>(
converted->getConvertedMatrix(),
w));
452template <
class Operator,
class Comm>
458template <
class Operator,
class Comm>
459PreconditionerFactory<Operator,Comm>&
460PreconditionerFactory<Operator,Comm>::instance()
462 static PreconditionerFactory singleton;
466template <
class Operator,
class Comm>
468PreconditionerFactory<Operator,Comm>::
469doCreate(
const Operator& op,
const PropertyTree& prm,
470 const std::function<Vector()> weightsCalculator,
471 std::size_t pressureIndex)
474 StandardPreconditioners<Operator,Comm>::add();
477 const std::string& type = prm.get<std::string>(
"type",
"ParOverILU0");
478 auto it = creators_.find(type);
479 if (it == creators_.end()) {
480 std::ostringstream msg;
481 msg <<
"Preconditioner type " << type <<
" is not registered in the factory. Available types are: ";
482 for (
const auto& prec : creators_) {
483 msg << prec.first <<
' ';
486 OPM_THROW(std::invalid_argument, msg.str());
488 return it->second(op, prm, weightsCalculator, pressureIndex);
491template <
class Operator,
class Comm>
493PreconditionerFactory<Operator,Comm>::
494doCreate(
const Operator& op,
const PropertyTree& prm,
495 const std::function<Vector()> weightsCalculator,
496 std::size_t pressureIndex,
const Comm& comm)
499 StandardPreconditioners<Operator,Comm>::add();
502 const std::string& type = prm.get<std::string>(
"type",
"ParOverILU0");
503 auto it = parallel_creators_.find(type);
504 if (it == parallel_creators_.end()) {
505 std::ostringstream msg;
506 msg <<
"Parallel preconditioner type " << type
507 <<
" is not registered in the factory. Available types are: ";
508 for (
const auto& prec : parallel_creators_) {
509 msg << prec.first <<
' ';
512 OPM_THROW(std::invalid_argument, msg.str());
514 return it->second(op, prm, weightsCalculator, pressureIndex, comm);
517template <
class Operator,
class Comm>
518void PreconditionerFactory<Operator,Comm>::
519doAddCreator(
const std::string& type, Creator c)
524template <
class Operator,
class Comm>
525void PreconditionerFactory<Operator,Comm>::
526doAddCreator(
const std::string& type, ParCreator c)
528 parallel_creators_[type] = c;
531template <
class Operator,
class Comm>
536 std::size_t pressureIndex)
541template <
class Operator,
class Comm>
546 std::size_t pressureIndex)
552template <
class Operator,
class Comm>
556 std::size_t pressureIndex)
558 return instance().doCreate(
op, prm, std::function<Vector()>(), pressureIndex, comm);
561template <
class Operator,
class Comm>
565 instance().doAddCreator(type,
creator);
568template <
class Operator,
class Comm>
572 instance().doAddCreator(type,
creator);
575using CommSeq = Dune::Amg::SequentialInformation;
578using OpFSeq = Dune::MatrixAdapter<Dune::BCRSMatrix<Dune::FieldMatrix<double,Dim,Dim>>,
579 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
580 Dune::BlockVector<Dune::FieldVector<double,Dim>>>;
582using OpBSeq = Dune::MatrixAdapter<Dune::BCRSMatrix<Opm::MatrixBlock<double,Dim,Dim>>,
583 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
584 Dune::BlockVector<Dune::FieldVector<double,Dim>>>;
586template<
int Dim,
bool overlap>
588 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
589 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
592template<
int Dim,
bool overlap>
594 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
595 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
599using CommPar = Dune::OwnerOverlapCopyCommunication<int,int>;
602using OpFPar = Dune::OverlappingSchwarzOperator<Dune::BCRSMatrix<Dune::FieldMatrix<double,Dim,Dim>>,
603 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
604 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
608using OpBPar = Dune::OverlappingSchwarzOperator<Dune::BCRSMatrix<MatrixBlock<double,Dim,Dim>>,
609 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
610 Dune::BlockVector<Dune::FieldVector<double,Dim>>,
613#define INSTANCE_PF_PAR(Dim) \
614 template class PreconditionerFactory<OpBSeq<Dim>,CommPar>; \
615 template class PreconditionerFactory<OpFPar<Dim>,CommPar>; \
616 template class PreconditionerFactory<OpBPar<Dim>,CommPar>; \
617 template class PreconditionerFactory<OpW<Dim,false>,CommPar>; \
618 template class PreconditionerFactory<OpWG<Dim,true>,CommPar>; \
619 template class PreconditionerFactory<OpBPar<Dim>,CommSeq>;
622#define INSTANCE_PF_SEQ(Dim) \
623 template class PreconditionerFactory<OpFSeq<Dim>,CommSeq>; \
624 template class PreconditionerFactory<OpBSeq<Dim>,CommSeq>; \
625 template class PreconditionerFactory<OpW<Dim,false>,CommSeq>; \
626 template class PreconditionerFactory<OpWG<Dim,true>,CommSeq>;
629#define INSTANCE_PF(Dim) \
630 INSTANCE_PF_PAR(Dim) \
633#define INSTANCE_PF(Dim) \
Parallel algebraic multigrid based on agglomeration.
Definition amgcpr.hh:88
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:534
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:563
std::shared_ptr< Dune::PreconditionerWithUpdate< Vector, Vector > > PrecPtr
The type of pointer returned by create().
Definition PreconditionerFactory.hpp:71
Definition PropertyTree.hpp:37
Sequential ILU0 preconditioner on the GPU through the CuSparse library.
Definition CuSeqILU0.hpp:50
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 BlackoilPhases.hpp:27
MILU_VARIANT
Definition MILU.hpp:34
@ ILU
Do not perform modified ILU.
Definition PreconditionerFactory.hpp:43
Definition PreconditionerFactory_impl.hpp:58
Definition PreconditionerFactory_impl.hpp:141
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:271