153 using namespace Dune;
157 using M =
typename F::Matrix;
158 using V =
typename F::Vector;
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);
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));
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));
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);
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());
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);
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);
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);
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);
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");
212 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
216 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
218 }
else if (smoother ==
"DILU") {
220 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
221 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
224 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
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;
232 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
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;
240 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
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;
248 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
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;
256 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
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;
264 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
267 OPM_THROW(std::invalid_argument,
"Properties: No smoother with name " + smoother +
".");
275 const std::function<V()> weightsCalculator,
276 std::size_t pressureIndex,
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");
283 using Scalar =
typename V::field_type;
285 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
286 op, prm, weightsCalculator, pressureIndex, comm);
288 F::addCreator(
"cprt",
291 const std::function<V()> weightsCalculator,
292 std::size_t pressureIndex,
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");
299 using Scalar =
typename V::field_type;
301 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
302 op, prm, weightsCalculator, pressureIndex, comm);
305 if constexpr (std::is_same_v<O, WellModelGhostLastMatrixAdapter<M, V, V, true>>) {
306 F::addCreator(
"cprw",
309 const std::function<V()> weightsCalculator,
310 std::size_t pressureIndex,
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");
317 using Scalar =
typename V::field_type;
319 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
320 op, prm, weightsCalculator, pressureIndex, comm);
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);
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);
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;
342 auto gpuJac = std::make_shared<GpuJac>(op.getmat(), w);
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);
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;
354 auto gpuDILU = std::make_shared<GpuDILU>(op.getmat(), split_matrix, tune_gpu_kernels);
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);
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;
367 auto gpuilu0 = std::make_shared<OpmGpuILU0>(op.getmat(), split_matrix, tune_gpu_kernels, store_factorization_as_float);
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);
378 createParILU(
const Operator& op,
const PropertyTree& prm,
const Comm& comm,
const int ilulevel)
381 using M =
typename F::Matrix;
382 using V =
typename F::Vector;
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);
390 return std::make_shared<ParallelOverlappingILU0<M, V, V, Comm>>(
391 op.getmat(), comm, w,
MILU_VARIANT::ILU, num_interior, redblack, reorder_spheres);
393 return std::make_shared<ParallelOverlappingILU0<M, V, V, Comm>>(
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())) {
410 highest_interior_index = std::max(highest_interior_index, ind.local().local());
413 if (highest_interior_index + 1 == interior_count) {
414 return interior_count;
425 using namespace Dune;
427 using C = Dune::Amg::SequentialInformation;
429 using M =
typename F::Matrix;
430 using V =
typename F::Vector;
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>>(
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);
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>>(
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>>(
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());
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);
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);
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);
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);
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>;
488 }
else if (smoother ==
"Jac") {
489 using Smoother = SeqJac<M, V, V>;
491 }
else if (smoother ==
"GS") {
492 using Smoother = SeqGS<M, V, V>;
494 }
else if (smoother ==
"DILU") {
497 }
else if (smoother ==
"SOR") {
498 using Smoother = SeqSOR<M, V, V>;
500 }
else if (smoother ==
"SSOR") {
501 using Smoother = SeqSSOR<M, V, V>;
503 }
else if (smoother ==
"ILUn") {
504 using Smoother = SeqILU<M, V, V>;
507 OPM_THROW(std::invalid_argument,
"Properties: No smoother with name " + smoother +
".");
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>;
515 }
else if (smoother ==
"Jac") {
516 using Smoother = SeqJac<M, V, V>;
518 }
else if (smoother ==
"SOR") {
519 using Smoother = SeqSOR<M, V, V>;
521 }
else if (smoother ==
"GS") {
522 using Smoother = SeqGS<M, V, V>;
524 }
else if (smoother ==
"SSOR") {
525 using Smoother = SeqSSOR<M, V, V>;
527 }
else if (smoother ==
"ILUn") {
528 using Smoother = SeqILU<M, V, V>;
531 OPM_THROW(std::invalid_argument,
"Properties: No smoother with name " + smoother +
".");
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");
540 Dune::Amg::Parameters parms;
541 parms.setNoPreSmoothSteps(1);
542 parms.setNoPostSmoothSteps(1);
543 return getRebuildOnUpdateWrapper<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
547 if constexpr (std::is_same_v<O, WellModelMatrixAdapter<M, V, V, false>>) {
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");
554 using Scalar =
typename V::field_type;
555 using LevelTransferPolicy
557 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
558 op, prm, weightsCalculator, pressureIndex);
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");
568 using Scalar =
typename V::field_type;
570 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
571 op, prm, weightsCalculator, pressureIndex);
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");
579 using Scalar =
typename V::field_type;
581 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
582 op, prm, weightsCalculator, pressureIndex);
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));
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);
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;
616 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, GPUJac>>(
617 std::make_shared<GPUJac>(op.getmat(), w));
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);
625 using field_type =
typename V::field_type;
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));
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;
636 return std::make_shared<gpuistl::PreconditionerAdapter<V, V, GPUDILU>>(std::make_shared<GPUDILU>(op.getmat(), split_matrix, tune_gpu_kernels));
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);
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>>;
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);