159 using namespace Dune;
163 using M =
typename F::Matrix;
164 using V =
typename F::Vector;
166 F::addCreator(
"ILU0", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
167 return createParILU(op, prm, comm, 0);
169 F::addCreator(
"ParOverILU0",
170 [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
171 return createParILU(op, prm, comm, prm.get<
int>(
"ilulevel", 0));
173 F::addCreator(
"ILUn", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
174 return createParILU(op, prm, comm, prm.get<
int>(
"ilulevel", 0));
176 F::addCreator(
"DuneILU", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
177 const int n = prm.get<
int>(
"ilulevel", 0);
178 const double w = prm.get<
double>(
"relaxation", 1.0);
179 const bool resort = prm.get<
bool>(
"resort",
false);
180 return wrapBlockPreconditioner<RebuildOnUpdatePreconditioner<Dune::SeqILU<M, V, V>>>(
181 comm, op.getmat(), n, w, resort);
183 F::addCreator(
"DILU", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
184 DUNE_UNUSED_PARAMETER(prm);
185 return wrapBlockPreconditioner<MultithreadDILU<M, V, V>>(comm, op.getmat());
187 F::addCreator(
"Jac", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
188 const int n = prm.get<
int>(
"repeats", 1);
189 const double w = prm.get<
double>(
"relaxation", 1.0);
190 return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqJac<M, V, V>>>(comm, op.getmat(), n, w);
192 F::addCreator(
"GS", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
193 const int n = prm.get<
int>(
"repeats", 1);
194 const double w = prm.get<
double>(
"relaxation", 1.0);
195 return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqGS<M, V, V>>>(comm, op.getmat(), n, w);
197 F::addCreator(
"SOR", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
198 const int n = prm.get<
int>(
"repeats", 1);
199 const double w = prm.get<
double>(
"relaxation", 1.0);
200 return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSOR<M, V, V>>>(comm, op.getmat(), n, w);
202 F::addCreator(
"SSOR", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
203 const int n = prm.get<
int>(
"repeats", 1);
204 const double w = prm.get<
double>(
"relaxation", 1.0);
205 return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSSOR<M, V, V>>>(comm, op.getmat(), n, w);
212 if constexpr (std::is_same_v<O, Dune::OverlappingSchwarzOperator<M, V, V, C>>) {
213 F::addCreator(
"amg", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
214 using PrecPtr = std::shared_ptr<Dune::PreconditionerWithUpdate<V, V>>;
215 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
217 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
221 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
223 }
else if (smoother ==
"DILU") {
225 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
226 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
229 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
231 }
else if (smoother ==
"Jac") {
232 using SeqSmoother = SeqJac<M, V, V>;
233 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
234 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
237 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
239 }
else if (smoother ==
"GS") {
240 using SeqSmoother = SeqGS<M, V, V>;
241 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
242 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
245 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
247 }
else if (smoother ==
"SOR") {
248 using SeqSmoother = SeqSOR<M, V, V>;
249 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
250 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
253 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
255 }
else if (smoother ==
"SSOR") {
256 using SeqSmoother = SeqSSOR<M, V, V>;
257 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
258 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
261 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
263 }
else if (smoother ==
"ILUn") {
264 using SeqSmoother = SeqILU<M, V, V>;
265 using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
266 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
269 PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
272 OPM_THROW(std::invalid_argument,
"Properties: No smoother with name " + smoother +
".");
280 const std::function<V()> weightsCalculator,
281 std::size_t pressureIndex,
283 assert(weightsCalculator);
284 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
285 OPM_THROW(std::logic_error,
286 "Pressure index out of bounds. It needs to specified for CPR");
289 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
290 op, prm, weightsCalculator, pressureIndex, comm);
292 F::addCreator(
"cprt",
295 const std::function<V()> weightsCalculator,
296 std::size_t pressureIndex,
298 assert(weightsCalculator);
299 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
300 OPM_THROW(std::logic_error,
301 "Pressure index out of bounds. It needs to specified for CPR");
304 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
305 op, prm, weightsCalculator, pressureIndex, comm);
308 if constexpr (std::is_same_v<O, WellModelGhostLastMatrixAdapter<M, V, V, true>>) {
309 F::addCreator(
"cprw",
312 const std::function<V()> weightsCalculator,
313 std::size_t pressureIndex,
315 assert(weightsCalculator);
316 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
317 OPM_THROW(std::logic_error,
318 "Pressure index out of bounds. It needs to specified for CPR");
321 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
322 op, prm, weightsCalculator, pressureIndex, comm);
327 F::addCreator(
"CUILU0", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
328 const double w = prm.get<
double>(
"relaxation", 1.0);
329 using field_type =
typename V::field_type;
330 using CuILU0 =
typename Opm::cuistl::
331 CuSeqILU0<M, Opm::cuistl::CuVector<field_type>, Opm::cuistl::CuVector<field_type>>;
332 auto cuILU0 = std::make_shared<CuILU0>(op.getmat(), w);
334 auto adapted = std::make_shared<Opm::cuistl::PreconditionerAdapter<V, V, CuILU0>>(cuILU0);
335 auto wrapped = std::make_shared<Opm::cuistl::CuBlockPreconditioner<V, V, Comm>>(adapted, comm);
339 F::addCreator(
"CUJac", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
340 const double w = prm.get<
double>(
"relaxation", 1.0);
341 using field_type =
typename V::field_type;
344 auto cuJac = std::make_shared<CuJac>(op.getmat(), w);
346 auto adapted = std::make_shared<Opm::cuistl::PreconditionerAdapter<V, V, CuJac>>(cuJac);
347 auto wrapped = std::make_shared<Opm::cuistl::CuBlockPreconditioner<V, V, Comm>>(adapted, comm);
351 F::addCreator(
"CUDILU", [](
const O& op, [[maybe_unused]]
const P& prm,
const std::function<V()>&, std::size_t,
const C& comm) {
352 using field_type =
typename V::field_type;
354 auto cuDILU = std::make_shared<CuDILU>(op.getmat());
356 auto adapted = std::make_shared<Opm::cuistl::PreconditionerAdapter<V, V, CuDILU>>(cuDILU);
357 auto wrapped = std::make_shared<Opm::cuistl::CuBlockPreconditioner<V, V, Comm>>(adapted, comm);
365 createParILU(
const Operator& op,
const PropertyTree& prm,
const Comm& comm,
const int ilulevel)
368 using M =
typename F::Matrix;
369 using V =
typename F::Vector;
371 const double w = prm.get<
double>(
"relaxation", 1.0);
372 const bool redblack = prm.get<
bool>(
"redblack",
false);
373 const bool reorder_spheres = prm.get<
bool>(
"reorder_spheres",
false);
377 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, Comm>>(
380 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, Comm>>(
391 std::size_t interior_count = 0;
392 std::size_t highest_interior_index = 0;
393 const auto& is = comm.indexSet();
394 for (
const auto& ind : is) {
395 if (Comm::OwnerSet::contains(ind.local().attribute())) {
397 highest_interior_index = std::max(highest_interior_index, ind.local().local());
400 if (highest_interior_index + 1 == interior_count) {
401 return interior_count;
412 using namespace Dune;
414 using C = Dune::Amg::SequentialInformation;
416 using M =
typename F::Matrix;
417 using V =
typename F::Vector;
419 F::addCreator(
"ILU0", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
420 const double w = prm.get<
double>(
"relaxation", 1.0);
421 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
424 F::addCreator(
"DuneILU", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
425 const double w = prm.get<
double>(
"relaxation", 1.0);
426 const int n = prm.get<
int>(
"ilulevel", 0);
427 const bool resort = prm.get<
bool>(
"resort",
false);
428 return getRebuildOnUpdateWrapper<Dune::SeqILU<M, V, V>>(op.getmat(), n, w, resort);
430 F::addCreator(
"ParOverILU0", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
431 const double w = prm.get<
double>(
"relaxation", 1.0);
432 const int n = prm.get<
int>(
"ilulevel", 0);
433 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
436 F::addCreator(
"ILUn", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
437 const int n = prm.get<
int>(
"ilulevel", 0);
438 const double w = prm.get<
double>(
"relaxation", 1.0);
439 return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
442 F::addCreator(
"DILU", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
443 DUNE_UNUSED_PARAMETER(prm);
444 return std::make_shared<MultithreadDILU<M, V, V>>(op.getmat());
446 F::addCreator(
"Jac", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
447 const int n = prm.get<
int>(
"repeats", 1);
448 const double w = prm.get<
double>(
"relaxation", 1.0);
449 return getDummyUpdateWrapper<SeqJac<M, V, V>>(op.getmat(), n, w);
451 F::addCreator(
"GS", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
452 const int n = prm.get<
int>(
"repeats", 1);
453 const double w = prm.get<
double>(
"relaxation", 1.0);
454 return getDummyUpdateWrapper<SeqGS<M, V, V>>(op.getmat(), n, w);
456 F::addCreator(
"SOR", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
457 const int n = prm.get<
int>(
"repeats", 1);
458 const double w = prm.get<
double>(
"relaxation", 1.0);
459 return getDummyUpdateWrapper<SeqSOR<M, V, V>>(op.getmat(), n, w);
461 F::addCreator(
"SSOR", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
462 const int n = prm.get<
int>(
"repeats", 1);
463 const double w = prm.get<
double>(
"relaxation", 1.0);
464 return getDummyUpdateWrapper<SeqSSOR<M, V, V>>(op.getmat(), n, w);
469 if constexpr (std::is_same_v<O, Dune::MatrixAdapter<M, V, V>>) {
470 F::addCreator(
"amg", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
471 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
472 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
473 using Smoother = SeqILU<M, V, V>;
475 }
else if (smoother ==
"Jac") {
476 using Smoother = SeqJac<M, V, V>;
478 }
else if (smoother ==
"GS") {
479 using Smoother = SeqGS<M, V, V>;
481 }
else if (smoother ==
"DILU") {
484 }
else if (smoother ==
"SOR") {
485 using Smoother = SeqSOR<M, V, V>;
487 }
else if (smoother ==
"SSOR") {
488 using Smoother = SeqSSOR<M, V, V>;
490 }
else if (smoother ==
"ILUn") {
491 using Smoother = SeqILU<M, V, V>;
494 OPM_THROW(std::invalid_argument,
"Properties: No smoother with name " + smoother +
".");
497 F::addCreator(
"kamg", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
498 const std::string smoother = prm.get<std::string>(
"smoother",
"ParOverILU0");
499 if (smoother ==
"ILU0" || smoother ==
"ParOverILU0") {
500 using Smoother = SeqILU<M, V, V>;
502 }
else if (smoother ==
"Jac") {
503 using Smoother = SeqJac<M, V, V>;
505 }
else if (smoother ==
"SOR") {
506 using Smoother = SeqSOR<M, V, V>;
508 }
else if (smoother ==
"GS") {
509 using Smoother = SeqGS<M, V, V>;
511 }
else if (smoother ==
"SSOR") {
512 using Smoother = SeqSSOR<M, V, V>;
514 }
else if (smoother ==
"ILUn") {
515 using Smoother = SeqILU<M, V, V>;
518 OPM_THROW(std::invalid_argument,
"Properties: No smoother with name " + smoother +
".");
521 F::addCreator(
"famg", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
523 Dune::Amg::Parameters parms;
524 parms.setNoPreSmoothSteps(1);
525 parms.setNoPostSmoothSteps(1);
526 return getRebuildOnUpdateWrapper<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
529 if constexpr (std::is_same_v<O, WellModelMatrixAdapter<M, V, V, false>>) {
532 [](
const O& op,
const P& prm,
const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
533 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
534 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
536 using LevelTransferPolicy
538 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
539 op, prm, weightsCalculator, pressureIndex);
545 [](
const O& op,
const P& prm,
const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
546 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
547 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
550 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
551 op, prm, weightsCalculator, pressureIndex);
555 [](
const O& op,
const P& prm,
const std::function<V()>& weightsCalculator, std::size_t pressureIndex) {
556 if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
557 OPM_THROW(std::logic_error,
"Pressure index out of bounds. It needs to specified for CPR");
560 return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(
561 op, prm, weightsCalculator, pressureIndex);
565 F::addCreator(
"CUILU0", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
566 const double w = prm.get<
double>(
"relaxation", 1.0);
567 using field_type =
typename V::field_type;
568 using CuILU0 =
typename Opm::cuistl::
569 CuSeqILU0<M, Opm::cuistl::CuVector<field_type>, Opm::cuistl::CuVector<field_type>>;
570 return std::make_shared<Opm::cuistl::PreconditionerAdapter<V, V, CuILU0>>(
571 std::make_shared<CuILU0>(op.getmat(), w));
574 F::addCreator(
"CUILU0Float", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
575 const double w = prm.get<
double>(
"relaxation", 1.0);
576 using block_type =
typename V::block_type;
577 using VTo = Dune::BlockVector<Dune::FieldVector<float, block_type::dimension>>;
578 using matrix_type_to =
579 typename Dune::BCRSMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
580 using CuILU0 =
typename Opm::cuistl::
581 CuSeqILU0<matrix_type_to, Opm::cuistl::CuVector<float>, Opm::cuistl::CuVector<float>>;
584 auto converted = std::make_shared<Converter>(op.getmat());
585 auto adapted = std::make_shared<Adapter>(std::make_shared<CuILU0>(converted->getConvertedMatrix(), w));
586 converted->setUnderlyingPreconditioner(adapted);
590 F::addCreator(
"CUJac", [](
const O& op,
const P& prm,
const std::function<V()>&, std::size_t) {
591 const double w = prm.get<
double>(
"relaxation", 1.0);
592 using field_type =
typename V::field_type;
595 return std::make_shared<Opm::cuistl::PreconditionerAdapter<V, V, CUJac>>(
596 std::make_shared<CUJac>(op.getmat(), w));
599 F::addCreator(
"CUDILU", [](
const O& op, [[maybe_unused]]
const P& prm,
const std::function<V()>&, std::size_t) {
600 using field_type =
typename V::field_type;
602 return std::make_shared<Opm::cuistl::PreconditionerAdapter<V, V, CUDILU>>(std::make_shared<CUDILU>(op.getmat()));
605 F::addCreator(
"CUDILUFloat", [](
const O& op, [[maybe_unused]]
const P& prm,
const std::function<V()>&, std::size_t) {
606 using block_type =
typename V::block_type;
607 using VTo = Dune::BlockVector<Dune::FieldVector<float, block_type::dimension>>;
608 using matrix_type_to =
typename Dune::BCRSMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
612 auto converted = std::make_shared<Converter>(op.getmat());
613 auto adapted = std::make_shared<Adapter>(std::make_shared<CuDILU>(converted->getConvertedMatrix()));
614 converted->setUnderlyingPreconditioner(adapted);