From d38aff6e25385c20bb7eb5c488a12e8353e2845b Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Tue, 29 Apr 2025 09:32:30 +0200 Subject: [PATCH] Splitting StandardPreconditioners to make it specializable for GPU --- CMakeLists.txt | 1 + CMakeLists_files.cmake | 7 + .../linalg/PreconditionerFactory_impl.hpp | 628 +----------------- .../linalg/StandardPreconditioners.hpp | 35 + .../StandardPreconditioners_gpu_mpi.hpp | 34 + .../StandardPreconditioners_gpu_serial.hpp | 35 + .../linalg/StandardPreconditioners_mpi.hpp | 399 +++++++++++ .../linalg/StandardPreconditioners_serial.hpp | 292 ++++++++ opm/simulators/linalg/is_gpu_operator.hpp | 55 ++ .../test_preconditioner_factory_gpu.cpp | 80 +++ 10 files changed, 940 insertions(+), 626 deletions(-) create mode 100644 opm/simulators/linalg/StandardPreconditioners.hpp create mode 100644 opm/simulators/linalg/StandardPreconditioners_gpu_mpi.hpp create mode 100644 opm/simulators/linalg/StandardPreconditioners_gpu_serial.hpp create mode 100644 opm/simulators/linalg/StandardPreconditioners_mpi.hpp create mode 100644 opm/simulators/linalg/StandardPreconditioners_serial.hpp create mode 100644 opm/simulators/linalg/is_gpu_operator.hpp create mode 100644 tests/gpuistl/test_preconditioner_factory_gpu.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index a4ec4e5cce..48a100d4e3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -756,6 +756,7 @@ if(CUDA_FOUND) gpu_smart_pointers is_gpu_pointer throw_macros_on_gpu + preconditioner_factory_gpu PROPERTIES LABELS ${gpu_label}) if(MPI_FOUND) set_tests_properties(GpuOwnerOverlapCopy PROPERTIES LABELS ${gpu_label}) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index a610098369..ba3d0c2f28 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -510,6 +510,8 @@ if (HAVE_CUDA) ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_gpu_resources.cu) ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_is_gpu_pointer.cpp) ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_throw_macros_on_gpu.cu) + ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_preconditioner_factory_gpu.cpp) + if(MPI_FOUND) ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_GpuOwnerOverlapCopy.cpp) endif() @@ -994,6 +996,11 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/PreconditionerFactoryGPUIncludeWrapper.hpp opm/simulators/linalg/PreconditionerFactory.hpp opm/simulators/linalg/PreconditionerFactory_impl.hpp + opm/simulators/linalg/StandardPreconditioners.hpp + opm/simulators/linalg/StandardPreconditioners_mpi.hpp + opm/simulators/linalg/StandardPreconditioners_serial.hpp + opm/simulators/linalg/StandardPreconditioners_gpu_serial.hpp + opm/simulators/linalg/StandardPreconditioners_gpu_mpi.hpp opm/simulators/linalg/PreconditionerWithUpdate.hpp opm/simulators/linalg/PressureBhpTransferPolicy.hpp opm/simulators/linalg/PressureSolverPolicy.hpp diff --git a/opm/simulators/linalg/PreconditionerFactory_impl.hpp b/opm/simulators/linalg/PreconditionerFactory_impl.hpp index 80fabb5d5b..c8ff504583 100644 --- a/opm/simulators/linalg/PreconditionerFactory_impl.hpp +++ b/opm/simulators/linalg/PreconditionerFactory_impl.hpp @@ -60,635 +60,11 @@ #include +#include + namespace Opm { -template -struct AMGSmootherArgsHelper -{ - static auto args(const PropertyTree& prm) - { - using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; - SmootherArgs smootherArgs; - smootherArgs.iterations = prm.get("iterations", 1); - // smootherArgs.overlap=SmootherArgs::vertex; - // smootherArgs.overlap=SmootherArgs::none; - // smootherArgs.overlap=SmootherArgs::aggregate; - smootherArgs.relaxationFactor = prm.get("relaxation", 1.0); - return smootherArgs; - } -}; -template -struct AMGSmootherArgsHelper> -{ - static auto args(const PropertyTree& prm) - { - using Smoother = ParallelOverlappingILU0; - using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; - SmootherArgs smootherArgs; - smootherArgs.iterations = prm.get("iterations", 1); - const int iluwitdh = prm.get("iluwidth", 0); - smootherArgs.setN(iluwitdh); - const MILU_VARIANT milu = convertString2Milu(prm.get("milutype", std::string("ilu"))); - smootherArgs.setMilu(milu); - // smootherArgs.overlap=SmootherArgs::vertex; - // smootherArgs.overlap=SmootherArgs::none; - // smootherArgs.overlap=SmootherArgs::aggregate; - smootherArgs.relaxationFactor = prm.get("relaxation", 1.0); - return smootherArgs; - } -}; - -// trailing return type with decltype used for detecting existence of setUseFixedOrder member function by overloading the setUseFixedOrder function -template -auto setUseFixedOrder(C& criterion, bool booleanValue) -> decltype(criterion.setUseFixedOrder(booleanValue)) -{ - return criterion.setUseFixedOrder(booleanValue); // Set flag to ensure that the matrices in the AMG hierarchy are constructed with deterministic indices. -} -template -void setUseFixedOrder(C&, ...) -{ - // do nothing, since the function setUseFixedOrder does not exist yet -} - -template -typename AMGHelper::Criterion -AMGHelper::criterion(const PropertyTree& prm) -{ - Criterion criterion(15, prm.get("coarsenTarget", 1200)); - criterion.setDefaultValuesIsotropic(2); - criterion.setAlpha(prm.get("alpha", 0.33)); - criterion.setBeta(prm.get("beta", 1e-5)); - criterion.setMaxLevel(prm.get("maxlevel", 15)); - criterion.setSkipIsolated(prm.get("skip_isolated", false)); - criterion.setNoPreSmoothSteps(prm.get("pre_smooth", 1)); - criterion.setNoPostSmoothSteps(prm.get("post_smooth", 1)); - criterion.setDebugLevel(prm.get("verbosity", 0)); - // As the default we request to accumulate data to 1 process always as our matrix - // graph might be unsymmetric and hence not supported by the PTScotch/ParMetis - // calls in DUNE. Accumulating to 1 skips PTScotch/ParMetis - criterion.setAccumulate(static_cast(prm.get("accumulate", 1))); - criterion.setProlongationDampingFactor(prm.get("prolongationdamping", 1.6)); - criterion.setMaxDistance(prm.get("maxdistance", 2)); - criterion.setMaxConnectivity(prm.get("maxconnectivity", 15)); - criterion.setMaxAggregateSize(prm.get("maxaggsize", 6)); - criterion.setMinAggregateSize(prm.get("minaggsize", 4)); - setUseFixedOrder(criterion, true); // If possible, set flag to ensure that the matrices in the AMG hierarchy are constructed with deterministic indices. - return criterion; -} - -template -template -typename AMGHelper::PrecPtr -AMGHelper::makeAmgPreconditioner(const Operator& op, - const PropertyTree& prm, - bool useKamg) -{ - auto crit = criterion(prm); - auto sargs = AMGSmootherArgsHelper::args(prm); - if (useKamg) { - using Type = Dune::DummyUpdatePreconditioner>; - return std::make_shared( - op, crit, sargs, prm.get("max_krylov", 1), prm.get("min_reduction", 1e-1)); - } else { - using Type = Dune::Amg::AMGCPR; - return std::make_shared(op, crit, sargs); - } -} - -template -struct StandardPreconditioners { - static void add() - { - using namespace Dune; - using O = Operator; - using C = Comm; - using F = PreconditionerFactory; - using M = typename F::Matrix; - using V = typename F::Vector; - using P = PropertyTree; - F::addCreator("ILU0", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - return createParILU(op, prm, comm, 0); - }); - F::addCreator("ParOverILU0", - [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - return createParILU(op, prm, comm, prm.get("ilulevel", 0)); - }); - F::addCreator("ILUn", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - return createParILU(op, prm, comm, prm.get("ilulevel", 0)); - }); - F::addCreator("DuneILU", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - const int n = prm.get("ilulevel", 0); - const double w = prm.get("relaxation", 1.0); - const bool resort = prm.get("resort", false); - return wrapBlockPreconditioner>>( - comm, op.getmat(), n, w, resort); - }); - F::addCreator("DILU", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - DUNE_UNUSED_PARAMETER(prm); - return wrapBlockPreconditioner>(comm, op.getmat()); - }); - F::addCreator("Jac", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - const int n = prm.get("repeats", 1); - const double w = prm.get("relaxation", 1.0); - return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); - }); - F::addCreator("GS", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - const int n = prm.get("repeats", 1); - const double w = prm.get("relaxation", 1.0); - return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); - }); - F::addCreator("SOR", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - const int n = prm.get("repeats", 1); - const double w = prm.get("relaxation", 1.0); - return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); - }); - F::addCreator("SSOR", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - const int n = prm.get("repeats", 1); - const double w = prm.get("relaxation", 1.0); - return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); - }); - - // Only add AMG preconditioners to the factory if the operator - // is the overlapping schwarz operator. This could be extended - // later, but at this point no other operators are compatible - // with the AMG hierarchy construction. - if constexpr (std::is_same_v> || - std::is_same_v>) { - F::addCreator("amg", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - using PrecPtr = std::shared_ptr>; - const std::string smoother = prm.get("smoother", "ParOverILU0"); - // TODO: merge this with ILUn, and possibly simplify the factory to only work with ILU? - if (smoother == "ILU0" || smoother == "ParOverILU0") { - using Smoother = ParallelOverlappingILU0; - auto crit = AMGHelper::criterion(prm); - auto sargs = AMGSmootherArgsHelper::args(prm); - PrecPtr prec = std::make_shared>(op, crit, sargs, comm); - return prec; - } else if (smoother == "DILU") { - using SeqSmoother = Dune::MultithreadDILU; - using Smoother = Dune::BlockPreconditioner; - using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; - SmootherArgs sargs; - auto crit = AMGHelper::criterion(prm); - PrecPtr prec = std::make_shared>(op, crit, sargs, comm); - return prec; - } else if (smoother == "Jac") { - using SeqSmoother = SeqJac; - using Smoother = Dune::BlockPreconditioner; - using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; - SmootherArgs sargs; - auto crit = AMGHelper::criterion(prm); - PrecPtr prec = std::make_shared>(op, crit, sargs, comm); - return prec; - } else if (smoother == "GS") { - using SeqSmoother = SeqGS; - using Smoother = Dune::BlockPreconditioner; - using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; - SmootherArgs sargs; - auto crit = AMGHelper::criterion(prm); - PrecPtr prec = std::make_shared>(op, crit, sargs, comm); - return prec; - } else if (smoother == "SOR") { - using SeqSmoother = SeqSOR; - using Smoother = Dune::BlockPreconditioner; - using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; - SmootherArgs sargs; - auto crit = AMGHelper::criterion(prm); - PrecPtr prec = std::make_shared>(op, crit, sargs, comm); - return prec; - } else if (smoother == "SSOR") { - using SeqSmoother = SeqSSOR; - using Smoother = Dune::BlockPreconditioner; - using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; - SmootherArgs sargs; - auto crit = AMGHelper::criterion(prm); - PrecPtr prec = std::make_shared>(op, crit, sargs, comm); - return prec; - } else if (smoother == "ILUn") { - using SeqSmoother = SeqILU; - using Smoother = Dune::BlockPreconditioner; - using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; - SmootherArgs sargs; - auto crit = AMGHelper::criterion(prm); - PrecPtr prec = std::make_shared>(op, crit, sargs, comm); - return prec; - } else { - OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + "."); - } - }); - } - - F::addCreator("cpr", - [](const O& op, - const P& prm, - const std::function weightsCalculator, - std::size_t pressureIndex, - const C& comm) { - assert(weightsCalculator); - if (pressureIndex == std::numeric_limits::max()) { - OPM_THROW(std::logic_error, - "Pressure index out of bounds. It needs to specified for CPR"); - } - using Scalar = typename V::field_type; - using LevelTransferPolicy = PressureTransferPolicy; - return std::make_shared>( - op, prm, weightsCalculator, pressureIndex, comm); - }); - F::addCreator("cprt", - [](const O& op, - const P& prm, - const std::function weightsCalculator, - std::size_t pressureIndex, - const C& comm) { - assert(weightsCalculator); - if (pressureIndex == std::numeric_limits::max()) { - OPM_THROW(std::logic_error, - "Pressure index out of bounds. It needs to specified for CPR"); - } - using Scalar = typename V::field_type; - using LevelTransferPolicy = PressureTransferPolicy; - return std::make_shared>( - op, prm, weightsCalculator, pressureIndex, comm); - }); - - if constexpr (std::is_same_v>) { - F::addCreator("cprw", - [](const O& op, - const P& prm, - const std::function weightsCalculator, - std::size_t pressureIndex, - const C& comm) { - assert(weightsCalculator); - if (pressureIndex == std::numeric_limits::max()) { - OPM_THROW(std::logic_error, - "Pressure index out of bounds. It needs to specified for CPR"); - } - using Scalar = typename V::field_type; - using LevelTransferPolicy = PressureBhpTransferPolicy; - return std::make_shared>( - op, prm, weightsCalculator, pressureIndex, comm); - }); - } - -#if HAVE_CUDA - F::addCreator("GPUILU0", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - const double w = prm.get("relaxation", 1.0); - using field_type = typename V::field_type; - using GpuILU0 = typename gpuistl:: - GpuSeqILU0, gpuistl::GpuVector>; - auto gpuILU0 = std::make_shared(op.getmat(), w); - - auto adapted = std::make_shared>(gpuILU0); - auto wrapped = std::make_shared>(adapted, comm); - return wrapped; - }); - - F::addCreator("GPUJAC", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { - const double w = prm.get("relaxation", 1.0); - using field_type = typename V::field_type; - using GpuJac = - typename gpuistl::GpuJac, gpuistl::GpuVector>; - auto gpuJac = std::make_shared(op.getmat(), w); - - auto adapted = std::make_shared>(gpuJac); - auto wrapped = std::make_shared>(adapted, comm); - return wrapped; - }); - - F::addCreator("GPUDILU", [](const O& op, [[maybe_unused]] const P& prm, const std::function&, std::size_t, const C& comm) { - const bool split_matrix = prm.get("split_matrix", true); - const bool tune_gpu_kernels = prm.get("tune_gpu_kernels", true); - const int mixed_precision_scheme = prm.get("mixed_precision_scheme", 0); - using field_type = typename V::field_type; - using GpuDILU = typename gpuistl::GpuDILU, gpuistl::GpuVector>; - auto gpuDILU = std::make_shared(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme); - - auto adapted = std::make_shared>(gpuDILU); - auto wrapped = std::make_shared>(adapted, comm); - return wrapped; - }); - - F::addCreator("OPMGPUILU0", [](const O& op, [[maybe_unused]] const P& prm, const std::function&, std::size_t, const C& comm) { - const bool split_matrix = prm.get("split_matrix", true); - const bool tune_gpu_kernels = prm.get("tune_gpu_kernels", true); - const int mixed_precision_scheme = prm.get("mixed_precision_scheme", 0); - using field_type = typename V::field_type; - using OpmGpuILU0 = typename gpuistl::OpmGpuILU0, gpuistl::GpuVector>; - auto gpuilu0 = std::make_shared(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme); - - auto adapted = std::make_shared>(gpuilu0); - auto wrapped = std::make_shared>(adapted, comm); - return wrapped; - }); -#endif - } - - - static typename PreconditionerFactory::PrecPtr - createParILU(const Operator& op, const PropertyTree& prm, const Comm& comm, const int ilulevel) - { - using F = PreconditionerFactory; - using M = typename F::Matrix; - using V = typename F::Vector; - - const double w = prm.get("relaxation", 1.0); - const bool redblack = prm.get("redblack", false); - const bool reorder_spheres = prm.get("reorder_spheres", false); - // Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner. - if (ilulevel == 0) { - const std::size_t num_interior = interiorIfGhostLast(comm); - assert(num_interior <= op.getmat().N()); - return std::make_shared>( - op.getmat(), comm, w, MILU_VARIANT::ILU, num_interior, redblack, reorder_spheres); - } else { - return std::make_shared>( - op.getmat(), comm, ilulevel, w, MILU_VARIANT::ILU, redblack, reorder_spheres); - } - } - - /// Helper method to determine if the local partitioning has the - /// K interior cells from [0, K-1] and ghost cells from [K, N-1]. - /// Returns K if true, otherwise returns N. This is motivated by - /// usage in the ParallelOverlappingILU0 preconditioner. - static std::size_t interiorIfGhostLast(const Comm& comm) - { - std::size_t interior_count = 0; - std::size_t highest_interior_index = 0; - const auto& is = comm.indexSet(); - for (const auto& ind : is) { - if (Comm::OwnerSet::contains(ind.local().attribute())) { - ++interior_count; - highest_interior_index = std::max(highest_interior_index, ind.local().local()); - } - } - if (highest_interior_index + 1 == interior_count) { - return interior_count; - } else { - return is.size(); - } - } -}; - -template -struct StandardPreconditioners { - static void add() - { - using namespace Dune; - using O = Operator; - using C = Dune::Amg::SequentialInformation; - using F = PreconditionerFactory; - using M = typename F::Matrix; - using V = typename F::Vector; - using P = PropertyTree; - F::addCreator("ILU0", [](const O& op, const P& prm, const std::function&, std::size_t) { - const double w = prm.get("relaxation", 1.0); - return std::make_shared>( - op.getmat(), 0, w, MILU_VARIANT::ILU); - }); - F::addCreator("DuneILU", [](const O& op, const P& prm, const std::function&, std::size_t) { - const double w = prm.get("relaxation", 1.0); - const int n = prm.get("ilulevel", 0); - const bool resort = prm.get("resort", false); - return getRebuildOnUpdateWrapper>(op.getmat(), n, w, resort); - }); - F::addCreator("ParOverILU0", [](const O& op, const P& prm, const std::function&, std::size_t) { - const double w = prm.get("relaxation", 1.0); - const int n = prm.get("ilulevel", 0); - return std::make_shared>( - op.getmat(), n, w, MILU_VARIANT::ILU); - }); - F::addCreator("ILUn", [](const O& op, const P& prm, const std::function&, std::size_t) { - const int n = prm.get("ilulevel", 0); - const double w = prm.get("relaxation", 1.0); - return std::make_shared>( - op.getmat(), n, w, MILU_VARIANT::ILU); - }); - F::addCreator("DILU", [](const O& op, const P& prm, const std::function&, std::size_t) { - DUNE_UNUSED_PARAMETER(prm); - return std::make_shared>(op.getmat()); - }); - F::addCreator("Jac", [](const O& op, const P& prm, const std::function&, std::size_t) { - const int n = prm.get("repeats", 1); - const double w = prm.get("relaxation", 1.0); - return getDummyUpdateWrapper>(op.getmat(), n, w); - }); - F::addCreator("GS", [](const O& op, const P& prm, const std::function&, std::size_t) { - const int n = prm.get("repeats", 1); - const double w = prm.get("relaxation", 1.0); - return getDummyUpdateWrapper>(op.getmat(), n, w); - }); - F::addCreator("SOR", [](const O& op, const P& prm, const std::function&, std::size_t) { - const int n = prm.get("repeats", 1); - const double w = prm.get("relaxation", 1.0); - return getDummyUpdateWrapper>(op.getmat(), n, w); - }); - F::addCreator("SSOR", [](const O& op, const P& prm, const std::function&, std::size_t) { - const int n = prm.get("repeats", 1); - const double w = prm.get("relaxation", 1.0); - return getDummyUpdateWrapper>(op.getmat(), n, w); - }); - - // Only add AMG preconditioners to the factory if the operator - // is an actual matrix operator. - if constexpr (std::is_same_v>) { - F::addCreator("amg", [](const O& op, const P& prm, const std::function&, std::size_t) { - const std::string smoother = prm.get("smoother", "ParOverILU0"); - if (smoother == "ILU0" || smoother == "ParOverILU0") { - using Smoother = SeqILU; - return AMGHelper::template makeAmgPreconditioner(op, prm); - } else if (smoother == "Jac") { - using Smoother = SeqJac; - return AMGHelper::template makeAmgPreconditioner(op, prm); - } else if (smoother == "GS") { - using Smoother = SeqGS; - return AMGHelper::template makeAmgPreconditioner(op, prm); - } else if (smoother == "DILU") { - using Smoother = MultithreadDILU; - return AMGHelper::template makeAmgPreconditioner(op, prm); - } else if (smoother == "SOR") { - using Smoother = SeqSOR; - return AMGHelper::template makeAmgPreconditioner(op, prm); - } else if (smoother == "SSOR") { - using Smoother = SeqSSOR; - return AMGHelper::template makeAmgPreconditioner(op, prm); - } else if (smoother == "ILUn") { - using Smoother = SeqILU; - return AMGHelper::template makeAmgPreconditioner(op, prm); - } else { - OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + "."); - } - }); - F::addCreator("kamg", [](const O& op, const P& prm, const std::function&, std::size_t) { - const std::string smoother = prm.get("smoother", "ParOverILU0"); - if (smoother == "ILU0" || smoother == "ParOverILU0") { - using Smoother = SeqILU; - return AMGHelper::template makeAmgPreconditioner(op, prm, true); - } else if (smoother == "Jac") { - using Smoother = SeqJac; - return AMGHelper::template makeAmgPreconditioner(op, prm, true); - } else if (smoother == "SOR") { - using Smoother = SeqSOR; - return AMGHelper::template makeAmgPreconditioner(op, prm, true); - } else if (smoother == "GS") { - using Smoother = SeqGS; - return AMGHelper::template makeAmgPreconditioner(op, prm, true); - } else if (smoother == "SSOR") { - using Smoother = SeqSSOR; - return AMGHelper::template makeAmgPreconditioner(op, prm, true); - } else if (smoother == "ILUn") { - using Smoother = SeqILU; - return AMGHelper::template makeAmgPreconditioner(op, prm, true); - } else { - OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + "."); - } - }); - F::addCreator("famg", [](const O& op, const P& prm, const std::function&, std::size_t) { - if constexpr (std::is_same_v) { - OPM_THROW(std::logic_error, "famg requires UMFPack which is not available for floats"); - return nullptr; - } else { - auto crit = AMGHelper::criterion(prm); - Dune::Amg::Parameters parms; - parms.setNoPreSmoothSteps(1); - parms.setNoPostSmoothSteps(1); - return getRebuildOnUpdateWrapper>(op, crit, parms); - } - }); - -#if HAVE_AMGX - // Only add AMGX for scalar matrices - if constexpr (M::block_type::rows == 1 && M::block_type::cols == 1) { - F::addCreator("amgx", [](const O& op, const P& prm, const std::function&, std::size_t) { - auto prm_copy = prm; - prm_copy.put("setup_frequency", Opm::Parameters::Get()); - return std::make_shared>(op.getmat(), prm_copy); - }); - } -#endif - -#if HAVE_HYPRE - // Only add Hypre for scalar matrices - if constexpr (M::block_type::rows == 1 && M::block_type::cols == 1 && - std::is_same_v) { - F::addCreator("hypre", [](const O& op, const P& prm, const std::function&, std::size_t) { - return std::make_shared>(op.getmat(), prm); - }); - } -#endif - } - if constexpr (std::is_same_v>) { - F::addCreator( - "cprw", - [](const O& op, const P& prm, const std::function& weightsCalculator, std::size_t pressureIndex) { - if (pressureIndex == std::numeric_limits::max()) { - OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR"); - } - using Scalar = typename V::field_type; - using LevelTransferPolicy - = PressureBhpTransferPolicy; - return std::make_shared>( - op, prm, weightsCalculator, pressureIndex); - }); - } - - F::addCreator( - "cpr", - [](const O& op, const P& prm, const std::function& weightsCalculator, std::size_t pressureIndex) { - if (pressureIndex == std::numeric_limits::max()) { - OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR"); - } - using Scalar = typename V::field_type; - using LevelTransferPolicy = PressureTransferPolicy; - return std::make_shared>( - op, prm, weightsCalculator, pressureIndex); - }); - F::addCreator( - "cprt", - [](const O& op, const P& prm, const std::function& weightsCalculator, std::size_t pressureIndex) { - if (pressureIndex == std::numeric_limits::max()) { - OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR"); - } - using Scalar = typename V::field_type; - using LevelTransferPolicy = PressureTransferPolicy; - return std::make_shared>( - op, prm, weightsCalculator, pressureIndex); - }); - -#if HAVE_CUDA - F::addCreator("GPUILU0", [](const O& op, const P& prm, const std::function&, std::size_t) { - const double w = prm.get("relaxation", 1.0); - using field_type = typename V::field_type; - using GpuILU0 = typename gpuistl:: - GpuSeqILU0, gpuistl::GpuVector>; - return std::make_shared>( - std::make_shared(op.getmat(), w)); - }); - - F::addCreator("GPUILU0Float", [](const O& op, const P& prm, const std::function&, std::size_t) { - const double w = prm.get("relaxation", 1.0); - using block_type = typename V::block_type; - using VTo = Dune::BlockVector>; - using matrix_type_to = - typename Dune::BCRSMatrix>; - using GpuILU0 = typename gpuistl:: - GpuSeqILU0, gpuistl::GpuVector>; - using Adapter = typename gpuistl::PreconditionerAdapter; - using Converter = typename gpuistl::PreconditionerConvertFieldTypeAdapter; - auto converted = std::make_shared(op.getmat()); - auto adapted = std::make_shared(std::make_shared(converted->getConvertedMatrix(), w)); - converted->setUnderlyingPreconditioner(adapted); - return converted; - }); - - F::addCreator("GPUJAC", [](const O& op, const P& prm, const std::function&, std::size_t) { - const double w = prm.get("relaxation", 1.0); - using field_type = typename V::field_type; - using GPUJac = - typename gpuistl::GpuJac, gpuistl::GpuVector>; - return std::make_shared>( - std::make_shared(op.getmat(), w)); - }); - - F::addCreator("OPMGPUILU0", [](const O& op, [[maybe_unused]] const P& prm, const std::function&, std::size_t) { - const bool split_matrix = prm.get("split_matrix", true); - const bool tune_gpu_kernels = prm.get("tune_gpu_kernels", true); - const int mixed_precision_scheme = prm.get("mixed_precision_scheme", 0); - - using field_type = typename V::field_type; - using GPUILU0 = typename gpuistl::OpmGpuILU0, gpuistl::GpuVector>; - - return std::make_shared>(std::make_shared(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme)); - }); - - F::addCreator("GPUDILU", [](const O& op, [[maybe_unused]] const P& prm, const std::function&, std::size_t) { - const bool split_matrix = prm.get("split_matrix", true); - const bool tune_gpu_kernels = prm.get("tune_gpu_kernels", true); - const int mixed_precision_scheme = prm.get("mixed_precision_scheme", 0); - using field_type = typename V::field_type; - using GPUDILU = typename gpuistl::GpuDILU, gpuistl::GpuVector>; - return std::make_shared>(std::make_shared(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme)); - }); - - F::addCreator("GPUDILUFloat", [](const O& op, [[maybe_unused]] const P& prm, const std::function&, std::size_t) { - const bool split_matrix = prm.get("split_matrix", true); - const bool tune_gpu_kernels = prm.get("tune_gpu_kernels", true); - const int mixed_precision_scheme = prm.get("mixed_precision_scheme", 0); - - using block_type = typename V::block_type; - using VTo = Dune::BlockVector>; - using matrix_type_to = typename Dune::BCRSMatrix>; - using GpuDILU = typename gpuistl::GpuDILU, gpuistl::GpuVector>; - using Adapter = typename gpuistl::PreconditionerAdapter; - using Converter = typename gpuistl::PreconditionerConvertFieldTypeAdapter; - auto converted = std::make_shared(op.getmat()); - auto adapted = std::make_shared(std::make_shared(converted->getConvertedMatrix(), split_matrix, tune_gpu_kernels, mixed_precision_scheme)); - converted->setUnderlyingPreconditioner(adapted); - return converted; - }); -#endif - } -}; template PreconditionerFactory::PreconditionerFactory() diff --git a/opm/simulators/linalg/StandardPreconditioners.hpp b/opm/simulators/linalg/StandardPreconditioners.hpp new file mode 100644 index 0000000000..0d9bf1e465 --- /dev/null +++ b/opm/simulators/linalg/StandardPreconditioners.hpp @@ -0,0 +1,35 @@ +/* + Copyright 2025 Equinor ASA + + This file is part of the Open Porous Media project (OPM). + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_STANDARDPRECONDITIONERS_HEADER +#define OPM_STANDARDPRECONDITIONERS_HEADER + +// This is a convenience header to include all standard preconditioners. +// It includes the serial and MPI versions of the standard preconditioners. + +// Note that you probably should not include this header directly, but rather +// include the PreconditionerFactory.hpp header, which will handle this. + +#include +#include +#include + +#if HAVE_CUDA +#include +#include +#endif + +#endif diff --git a/opm/simulators/linalg/StandardPreconditioners_gpu_mpi.hpp b/opm/simulators/linalg/StandardPreconditioners_gpu_mpi.hpp new file mode 100644 index 0000000000..b641d3d9c6 --- /dev/null +++ b/opm/simulators/linalg/StandardPreconditioners_gpu_mpi.hpp @@ -0,0 +1,34 @@ +/* + Copyright 2025 Equinor ASA + + This file is part of the Open Porous Media project (OPM). + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_STANDARDPRECONDITIONERS_GPU_MPI_HEADER +#define OPM_STANDARDPRECONDITIONERS_GPU_MPI_HEADER + + +namespace Opm { + + +template +struct StandardPreconditioners>> { + static void add() + { + // No standard preconditioners for this type of operator. + } +}; + + +} +#endif diff --git a/opm/simulators/linalg/StandardPreconditioners_gpu_serial.hpp b/opm/simulators/linalg/StandardPreconditioners_gpu_serial.hpp new file mode 100644 index 0000000000..bb7b0d5109 --- /dev/null +++ b/opm/simulators/linalg/StandardPreconditioners_gpu_serial.hpp @@ -0,0 +1,35 @@ +/* + Copyright 2025 Equinor ASA + + This file is part of the Open Porous Media project (OPM). + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_STANDARDPRECONDITIONERS_GPU_SERIAL_HEADER +#define OPM_STANDARDPRECONDITIONERS_GPU_SERIAL_HEADER + + +namespace Opm { + + +template +struct StandardPreconditioners>> { + static void add() + { + // No standard preconditioners for this type of operator. + } +}; + + +} + +#endif diff --git a/opm/simulators/linalg/StandardPreconditioners_mpi.hpp b/opm/simulators/linalg/StandardPreconditioners_mpi.hpp new file mode 100644 index 0000000000..33d2bf4e0a --- /dev/null +++ b/opm/simulators/linalg/StandardPreconditioners_mpi.hpp @@ -0,0 +1,399 @@ +/* + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2019 SINTEF Digital, Mathematics and Cybernetics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_STANDARDPRECONDITIONERS_MPI_HEADER +#define OPM_STANDARDPRECONDITIONERS_MPI_HEADER + + +namespace Opm { + + +template +struct AMGSmootherArgsHelper +{ + static auto args(const PropertyTree& prm) + { + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + SmootherArgs smootherArgs; + smootherArgs.iterations = prm.get("iterations", 1); + // smootherArgs.overlap=SmootherArgs::vertex; + // smootherArgs.overlap=SmootherArgs::none; + // smootherArgs.overlap=SmootherArgs::aggregate; + smootherArgs.relaxationFactor = prm.get("relaxation", 1.0); + return smootherArgs; + } +}; + +template +struct AMGSmootherArgsHelper> +{ + static auto args(const PropertyTree& prm) + { + using Smoother = ParallelOverlappingILU0; + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + SmootherArgs smootherArgs; + smootherArgs.iterations = prm.get("iterations", 1); + const int iluwitdh = prm.get("iluwidth", 0); + smootherArgs.setN(iluwitdh); + const MILU_VARIANT milu = convertString2Milu(prm.get("milutype", std::string("ilu"))); + smootherArgs.setMilu(milu); + // smootherArgs.overlap=SmootherArgs::vertex; + // smootherArgs.overlap=SmootherArgs::none; + // smootherArgs.overlap=SmootherArgs::aggregate; + smootherArgs.relaxationFactor = prm.get("relaxation", 1.0); + return smootherArgs; + } +}; + +// trailing return type with decltype used for detecting existence of setUseFixedOrder member function by overloading the setUseFixedOrder function +template +auto setUseFixedOrder(C& criterion, bool booleanValue) -> decltype(criterion.setUseFixedOrder(booleanValue)) +{ + return criterion.setUseFixedOrder(booleanValue); // Set flag to ensure that the matrices in the AMG hierarchy are constructed with deterministic indices. +} +template +void setUseFixedOrder(C&, ...) +{ + // do nothing, since the function setUseFixedOrder does not exist yet +} + +template +typename AMGHelper::Criterion +AMGHelper::criterion(const PropertyTree& prm) +{ + Criterion criterion(15, prm.get("coarsenTarget", 1200)); + criterion.setDefaultValuesIsotropic(2); + criterion.setAlpha(prm.get("alpha", 0.33)); + criterion.setBeta(prm.get("beta", 1e-5)); + criterion.setMaxLevel(prm.get("maxlevel", 15)); + criterion.setSkipIsolated(prm.get("skip_isolated", false)); + criterion.setNoPreSmoothSteps(prm.get("pre_smooth", 1)); + criterion.setNoPostSmoothSteps(prm.get("post_smooth", 1)); + criterion.setDebugLevel(prm.get("verbosity", 0)); + // As the default we request to accumulate data to 1 process always as our matrix + // graph might be unsymmetric and hence not supported by the PTScotch/ParMetis + // calls in DUNE. Accumulating to 1 skips PTScotch/ParMetis + criterion.setAccumulate(static_cast(prm.get("accumulate", 1))); + criterion.setProlongationDampingFactor(prm.get("prolongationdamping", 1.6)); + criterion.setMaxDistance(prm.get("maxdistance", 2)); + criterion.setMaxConnectivity(prm.get("maxconnectivity", 15)); + criterion.setMaxAggregateSize(prm.get("maxaggsize", 6)); + criterion.setMinAggregateSize(prm.get("minaggsize", 4)); + setUseFixedOrder(criterion, true); // If possible, set flag to ensure that the matrices in the AMG hierarchy are constructed with deterministic indices. + return criterion; +} + +template +template +typename AMGHelper::PrecPtr +AMGHelper::makeAmgPreconditioner(const Operator& op, + const PropertyTree& prm, + bool useKamg) +{ + auto crit = criterion(prm); + auto sargs = AMGSmootherArgsHelper::args(prm); + if (useKamg) { + using Type = Dune::DummyUpdatePreconditioner>; + return std::make_shared( + op, crit, sargs, prm.get("max_krylov", 1), prm.get("min_reduction", 1e-1)); + } else { + using Type = Dune::Amg::AMGCPR; + return std::make_shared(op, crit, sargs); + } +} + +template // Note: Last argument is to allow partial specialization for GPU +struct StandardPreconditioners { + static void add() + { + using namespace Dune; + using O = Operator; + using C = Comm; + using F = PreconditionerFactory; + using M = typename F::Matrix; + using V = typename F::Vector; + using P = PropertyTree; + F::addCreator("ILU0", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + return createParILU(op, prm, comm, 0); + }); + F::addCreator("ParOverILU0", + [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + return createParILU(op, prm, comm, prm.get("ilulevel", 0)); + }); + F::addCreator("ILUn", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + return createParILU(op, prm, comm, prm.get("ilulevel", 0)); + }); + F::addCreator("DuneILU", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + const int n = prm.get("ilulevel", 0); + const double w = prm.get("relaxation", 1.0); + const bool resort = prm.get("resort", false); + return wrapBlockPreconditioner>>( + comm, op.getmat(), n, w, resort); + }); + F::addCreator("DILU", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + DUNE_UNUSED_PARAMETER(prm); + return wrapBlockPreconditioner>(comm, op.getmat()); + }); + F::addCreator("Jac", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + const int n = prm.get("repeats", 1); + const double w = prm.get("relaxation", 1.0); + return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); + }); + F::addCreator("GS", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + const int n = prm.get("repeats", 1); + const double w = prm.get("relaxation", 1.0); + return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); + }); + F::addCreator("SOR", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + const int n = prm.get("repeats", 1); + const double w = prm.get("relaxation", 1.0); + return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); + }); + F::addCreator("SSOR", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + const int n = prm.get("repeats", 1); + const double w = prm.get("relaxation", 1.0); + return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); + }); + + // Only add AMG preconditioners to the factory if the operator + // is the overlapping schwarz operator. This could be extended + // later, but at this point no other operators are compatible + // with the AMG hierarchy construction. + if constexpr (std::is_same_v> || + std::is_same_v>) { + F::addCreator("amg", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + using PrecPtr = std::shared_ptr>; + const std::string smoother = prm.get("smoother", "ParOverILU0"); + // TODO: merge this with ILUn, and possibly simplify the factory to only work with ILU? + if (smoother == "ILU0" || smoother == "ParOverILU0") { + using Smoother = ParallelOverlappingILU0; + auto crit = AMGHelper::criterion(prm); + auto sargs = AMGSmootherArgsHelper::args(prm); + PrecPtr prec = std::make_shared>(op, crit, sargs, comm); + return prec; + } else if (smoother == "DILU") { + using SeqSmoother = Dune::MultithreadDILU; + using Smoother = Dune::BlockPreconditioner; + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + SmootherArgs sargs; + auto crit = AMGHelper::criterion(prm); + PrecPtr prec = std::make_shared>(op, crit, sargs, comm); + return prec; + } else if (smoother == "Jac") { + using SeqSmoother = SeqJac; + using Smoother = Dune::BlockPreconditioner; + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + SmootherArgs sargs; + auto crit = AMGHelper::criterion(prm); + PrecPtr prec = std::make_shared>(op, crit, sargs, comm); + return prec; + } else if (smoother == "GS") { + using SeqSmoother = SeqGS; + using Smoother = Dune::BlockPreconditioner; + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + SmootherArgs sargs; + auto crit = AMGHelper::criterion(prm); + PrecPtr prec = std::make_shared>(op, crit, sargs, comm); + return prec; + } else if (smoother == "SOR") { + using SeqSmoother = SeqSOR; + using Smoother = Dune::BlockPreconditioner; + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + SmootherArgs sargs; + auto crit = AMGHelper::criterion(prm); + PrecPtr prec = std::make_shared>(op, crit, sargs, comm); + return prec; + } else if (smoother == "SSOR") { + using SeqSmoother = SeqSSOR; + using Smoother = Dune::BlockPreconditioner; + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + SmootherArgs sargs; + auto crit = AMGHelper::criterion(prm); + PrecPtr prec = std::make_shared>(op, crit, sargs, comm); + return prec; + } else if (smoother == "ILUn") { + using SeqSmoother = SeqILU; + using Smoother = Dune::BlockPreconditioner; + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + SmootherArgs sargs; + auto crit = AMGHelper::criterion(prm); + PrecPtr prec = std::make_shared>(op, crit, sargs, comm); + return prec; + } else { + OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + "."); + } + }); + } + + F::addCreator("cpr", + [](const O& op, + const P& prm, + const std::function weightsCalculator, + std::size_t pressureIndex, + const C& comm) { + assert(weightsCalculator); + if (pressureIndex == std::numeric_limits::max()) { + OPM_THROW(std::logic_error, + "Pressure index out of bounds. It needs to specified for CPR"); + } + using Scalar = typename V::field_type; + using LevelTransferPolicy = PressureTransferPolicy; + return std::make_shared>( + op, prm, weightsCalculator, pressureIndex, comm); + }); + F::addCreator("cprt", + [](const O& op, + const P& prm, + const std::function weightsCalculator, + std::size_t pressureIndex, + const C& comm) { + assert(weightsCalculator); + if (pressureIndex == std::numeric_limits::max()) { + OPM_THROW(std::logic_error, + "Pressure index out of bounds. It needs to specified for CPR"); + } + using Scalar = typename V::field_type; + using LevelTransferPolicy = PressureTransferPolicy; + return std::make_shared>( + op, prm, weightsCalculator, pressureIndex, comm); + }); + + if constexpr (std::is_same_v>) { + F::addCreator("cprw", + [](const O& op, + const P& prm, + const std::function weightsCalculator, + std::size_t pressureIndex, + const C& comm) { + assert(weightsCalculator); + if (pressureIndex == std::numeric_limits::max()) { + OPM_THROW(std::logic_error, + "Pressure index out of bounds. It needs to specified for CPR"); + } + using Scalar = typename V::field_type; + using LevelTransferPolicy = PressureBhpTransferPolicy; + return std::make_shared>( + op, prm, weightsCalculator, pressureIndex, comm); + }); + } + +#if HAVE_CUDA + F::addCreator("GPUILU0", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + const double w = prm.get("relaxation", 1.0); + using field_type = typename V::field_type; + using GpuILU0 = typename gpuistl:: + GpuSeqILU0, gpuistl::GpuVector>; + auto gpuILU0 = std::make_shared(op.getmat(), w); + + auto adapted = std::make_shared>(gpuILU0); + auto wrapped = std::make_shared>(adapted, comm); + return wrapped; + }); + + F::addCreator("GPUJAC", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { + const double w = prm.get("relaxation", 1.0); + using field_type = typename V::field_type; + using GpuJac = + typename gpuistl::GpuJac, gpuistl::GpuVector>; + auto gpuJac = std::make_shared(op.getmat(), w); + + auto adapted = std::make_shared>(gpuJac); + auto wrapped = std::make_shared>(adapted, comm); + return wrapped; + }); + + F::addCreator("GPUDILU", [](const O& op, [[maybe_unused]] const P& prm, const std::function&, std::size_t, const C& comm) { + const bool split_matrix = prm.get("split_matrix", true); + const bool tune_gpu_kernels = prm.get("tune_gpu_kernels", true); + const int mixed_precision_scheme = prm.get("mixed_precision_scheme", 0); + using field_type = typename V::field_type; + using GpuDILU = typename gpuistl::GpuDILU, gpuistl::GpuVector>; + auto gpuDILU = std::make_shared(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme); + + auto adapted = std::make_shared>(gpuDILU); + auto wrapped = std::make_shared>(adapted, comm); + return wrapped; + }); + + F::addCreator("OPMGPUILU0", [](const O& op, [[maybe_unused]] const P& prm, const std::function&, std::size_t, const C& comm) { + const bool split_matrix = prm.get("split_matrix", true); + const bool tune_gpu_kernels = prm.get("tune_gpu_kernels", true); + const int mixed_precision_scheme = prm.get("mixed_precision_scheme", 0); + using field_type = typename V::field_type; + using OpmGpuILU0 = typename gpuistl::OpmGpuILU0, gpuistl::GpuVector>; + auto gpuilu0 = std::make_shared(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme); + + auto adapted = std::make_shared>(gpuilu0); + auto wrapped = std::make_shared>(adapted, comm); + return wrapped; + }); +#endif + } + + + static typename PreconditionerFactory::PrecPtr + createParILU(const Operator& op, const PropertyTree& prm, const Comm& comm, const int ilulevel) + { + using F = PreconditionerFactory; + using M = typename F::Matrix; + using V = typename F::Vector; + + const double w = prm.get("relaxation", 1.0); + const bool redblack = prm.get("redblack", false); + const bool reorder_spheres = prm.get("reorder_spheres", false); + // Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner. + if (ilulevel == 0) { + const std::size_t num_interior = interiorIfGhostLast(comm); + assert(num_interior <= op.getmat().N()); + return std::make_shared>( + op.getmat(), comm, w, MILU_VARIANT::ILU, num_interior, redblack, reorder_spheres); + } else { + return std::make_shared>( + op.getmat(), comm, ilulevel, w, MILU_VARIANT::ILU, redblack, reorder_spheres); + } + } + + /// Helper method to determine if the local partitioning has the + /// K interior cells from [0, K-1] and ghost cells from [K, N-1]. + /// Returns K if true, otherwise returns N. This is motivated by + /// usage in the ParallelOverlappingILU0 preconditioner. + static std::size_t interiorIfGhostLast(const Comm& comm) + { + std::size_t interior_count = 0; + std::size_t highest_interior_index = 0; + const auto& is = comm.indexSet(); + for (const auto& ind : is) { + if (Comm::OwnerSet::contains(ind.local().attribute())) { + ++interior_count; + highest_interior_index = std::max(highest_interior_index, ind.local().local()); + } + } + if (highest_interior_index + 1 == interior_count) { + return interior_count; + } else { + return is.size(); + } + } +}; + + +} + +#endif \ No newline at end of file diff --git a/opm/simulators/linalg/StandardPreconditioners_serial.hpp b/opm/simulators/linalg/StandardPreconditioners_serial.hpp new file mode 100644 index 0000000000..20e83dd293 --- /dev/null +++ b/opm/simulators/linalg/StandardPreconditioners_serial.hpp @@ -0,0 +1,292 @@ +/* + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2019 SINTEF Digital, Mathematics and Cybernetics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_STANDARDPRECONDITIONERS_SERIAL_HPP +#define OPM_STANDARDPRECONDITIONERS_SERIAL_HPP + + +namespace Opm { + + +template +struct StandardPreconditioners>> { + static void add() + { + using namespace Dune; + using O = Operator; + using C = Dune::Amg::SequentialInformation; + using F = PreconditionerFactory; + using M = typename F::Matrix; + using V = typename F::Vector; + using P = PropertyTree; + F::addCreator("ILU0", [](const O& op, const P& prm, const std::function&, std::size_t) { + const double w = prm.get("relaxation", 1.0); + return std::make_shared>( + op.getmat(), 0, w, MILU_VARIANT::ILU); + }); + F::addCreator("DuneILU", [](const O& op, const P& prm, const std::function&, std::size_t) { + const double w = prm.get("relaxation", 1.0); + const int n = prm.get("ilulevel", 0); + const bool resort = prm.get("resort", false); + return getRebuildOnUpdateWrapper>(op.getmat(), n, w, resort); + }); + F::addCreator("ParOverILU0", [](const O& op, const P& prm, const std::function&, std::size_t) { + const double w = prm.get("relaxation", 1.0); + const int n = prm.get("ilulevel", 0); + return std::make_shared>( + op.getmat(), n, w, MILU_VARIANT::ILU); + }); + F::addCreator("ILUn", [](const O& op, const P& prm, const std::function&, std::size_t) { + const int n = prm.get("ilulevel", 0); + const double w = prm.get("relaxation", 1.0); + return std::make_shared>( + op.getmat(), n, w, MILU_VARIANT::ILU); + }); + F::addCreator("DILU", [](const O& op, const P& prm, const std::function&, std::size_t) { + DUNE_UNUSED_PARAMETER(prm); + return std::make_shared>(op.getmat()); + }); + F::addCreator("Jac", [](const O& op, const P& prm, const std::function&, std::size_t) { + const int n = prm.get("repeats", 1); + const double w = prm.get("relaxation", 1.0); + return getDummyUpdateWrapper>(op.getmat(), n, w); + }); + F::addCreator("GS", [](const O& op, const P& prm, const std::function&, std::size_t) { + const int n = prm.get("repeats", 1); + const double w = prm.get("relaxation", 1.0); + return getDummyUpdateWrapper>(op.getmat(), n, w); + }); + F::addCreator("SOR", [](const O& op, const P& prm, const std::function&, std::size_t) { + const int n = prm.get("repeats", 1); + const double w = prm.get("relaxation", 1.0); + return getDummyUpdateWrapper>(op.getmat(), n, w); + }); + F::addCreator("SSOR", [](const O& op, const P& prm, const std::function&, std::size_t) { + const int n = prm.get("repeats", 1); + const double w = prm.get("relaxation", 1.0); + return getDummyUpdateWrapper>(op.getmat(), n, w); + }); + + // Only add AMG preconditioners to the factory if the operator + // is an actual matrix operator. + if constexpr (std::is_same_v>) { + F::addCreator("amg", [](const O& op, const P& prm, const std::function&, std::size_t) { + const std::string smoother = prm.get("smoother", "ParOverILU0"); + if (smoother == "ILU0" || smoother == "ParOverILU0") { + using Smoother = SeqILU; + return AMGHelper::template makeAmgPreconditioner(op, prm); + } else if (smoother == "Jac") { + using Smoother = SeqJac; + return AMGHelper::template makeAmgPreconditioner(op, prm); + } else if (smoother == "GS") { + using Smoother = SeqGS; + return AMGHelper::template makeAmgPreconditioner(op, prm); + } else if (smoother == "DILU") { + using Smoother = MultithreadDILU; + return AMGHelper::template makeAmgPreconditioner(op, prm); + } else if (smoother == "SOR") { + using Smoother = SeqSOR; + return AMGHelper::template makeAmgPreconditioner(op, prm); + } else if (smoother == "SSOR") { + using Smoother = SeqSSOR; + return AMGHelper::template makeAmgPreconditioner(op, prm); + } else if (smoother == "ILUn") { + using Smoother = SeqILU; + return AMGHelper::template makeAmgPreconditioner(op, prm); + } else { + OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + "."); + } + }); + F::addCreator("kamg", [](const O& op, const P& prm, const std::function&, std::size_t) { + const std::string smoother = prm.get("smoother", "ParOverILU0"); + if (smoother == "ILU0" || smoother == "ParOverILU0") { + using Smoother = SeqILU; + return AMGHelper::template makeAmgPreconditioner(op, prm, true); + } else if (smoother == "Jac") { + using Smoother = SeqJac; + return AMGHelper::template makeAmgPreconditioner(op, prm, true); + } else if (smoother == "SOR") { + using Smoother = SeqSOR; + return AMGHelper::template makeAmgPreconditioner(op, prm, true); + } else if (smoother == "GS") { + using Smoother = SeqGS; + return AMGHelper::template makeAmgPreconditioner(op, prm, true); + } else if (smoother == "SSOR") { + using Smoother = SeqSSOR; + return AMGHelper::template makeAmgPreconditioner(op, prm, true); + } else if (smoother == "ILUn") { + using Smoother = SeqILU; + return AMGHelper::template makeAmgPreconditioner(op, prm, true); + } else { + OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + "."); + } + }); + F::addCreator("famg", [](const O& op, const P& prm, const std::function&, std::size_t) { + if constexpr (std::is_same_v) { + OPM_THROW(std::logic_error, "famg requires UMFPack which is not available for floats"); + return nullptr; + } else { + auto crit = AMGHelper::criterion(prm); + Dune::Amg::Parameters parms; + parms.setNoPreSmoothSteps(1); + parms.setNoPostSmoothSteps(1); + return getRebuildOnUpdateWrapper>(op, crit, parms); + } + }); + +#if HAVE_AMGX + // Only add AMGX for scalar matrices + if constexpr (M::block_type::rows == 1 && M::block_type::cols == 1) { + F::addCreator("amgx", [](const O& op, const P& prm, const std::function&, std::size_t) { + auto prm_copy = prm; + prm_copy.put("setup_frequency", Opm::Parameters::Get()); + return std::make_shared>(op.getmat(), prm_copy); + }); + } +#endif + +#if HAVE_HYPRE + // Only add Hypre for scalar matrices + if constexpr (M::block_type::rows == 1 && M::block_type::cols == 1 && + std::is_same_v) { + F::addCreator("hypre", [](const O& op, const P& prm, const std::function&, std::size_t) { + return std::make_shared>(op.getmat(), prm); + }); + } +#endif + } + if constexpr (std::is_same_v>) { + F::addCreator( + "cprw", + [](const O& op, const P& prm, const std::function& weightsCalculator, std::size_t pressureIndex) { + if (pressureIndex == std::numeric_limits::max()) { + OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR"); + } + using Scalar = typename V::field_type; + using LevelTransferPolicy + = PressureBhpTransferPolicy; + return std::make_shared>( + op, prm, weightsCalculator, pressureIndex); + }); + } + + F::addCreator( + "cpr", + [](const O& op, const P& prm, const std::function& weightsCalculator, std::size_t pressureIndex) { + if (pressureIndex == std::numeric_limits::max()) { + OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR"); + } + using Scalar = typename V::field_type; + using LevelTransferPolicy = PressureTransferPolicy; + return std::make_shared>( + op, prm, weightsCalculator, pressureIndex); + }); + F::addCreator( + "cprt", + [](const O& op, const P& prm, const std::function& weightsCalculator, std::size_t pressureIndex) { + if (pressureIndex == std::numeric_limits::max()) { + OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR"); + } + using Scalar = typename V::field_type; + using LevelTransferPolicy = PressureTransferPolicy; + return std::make_shared>( + op, prm, weightsCalculator, pressureIndex); + }); + +#if HAVE_CUDA + F::addCreator("GPUILU0", [](const O& op, const P& prm, const std::function&, std::size_t) { + const double w = prm.get("relaxation", 1.0); + using field_type = typename V::field_type; + using GpuILU0 = typename gpuistl:: + GpuSeqILU0, gpuistl::GpuVector>; + return std::make_shared>( + std::make_shared(op.getmat(), w)); + }); + + F::addCreator("GPUILU0Float", [](const O& op, const P& prm, const std::function&, std::size_t) { + const double w = prm.get("relaxation", 1.0); + using block_type = typename V::block_type; + using VTo = Dune::BlockVector>; + using matrix_type_to = + typename Dune::BCRSMatrix>; + using GpuILU0 = typename gpuistl:: + GpuSeqILU0, gpuistl::GpuVector>; + using Adapter = typename gpuistl::PreconditionerAdapter; + using Converter = typename gpuistl::PreconditionerConvertFieldTypeAdapter; + auto converted = std::make_shared(op.getmat()); + auto adapted = std::make_shared(std::make_shared(converted->getConvertedMatrix(), w)); + converted->setUnderlyingPreconditioner(adapted); + return converted; + }); + + F::addCreator("GPUJAC", [](const O& op, const P& prm, const std::function&, std::size_t) { + const double w = prm.get("relaxation", 1.0); + using field_type = typename V::field_type; + using GPUJac = + typename gpuistl::GpuJac, gpuistl::GpuVector>; + return std::make_shared>( + std::make_shared(op.getmat(), w)); + }); + + F::addCreator("OPMGPUILU0", [](const O& op, [[maybe_unused]] const P& prm, const std::function&, std::size_t) { + const bool split_matrix = prm.get("split_matrix", true); + const bool tune_gpu_kernels = prm.get("tune_gpu_kernels", true); + const int mixed_precision_scheme = prm.get("mixed_precision_scheme", 0); + + using field_type = typename V::field_type; + using GPUILU0 = typename gpuistl::OpmGpuILU0, gpuistl::GpuVector>; + + return std::make_shared>(std::make_shared(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme)); + }); + + F::addCreator("GPUDILU", [](const O& op, [[maybe_unused]] const P& prm, const std::function&, std::size_t) { + const bool split_matrix = prm.get("split_matrix", true); + const bool tune_gpu_kernels = prm.get("tune_gpu_kernels", true); + const int mixed_precision_scheme = prm.get("mixed_precision_scheme", 0); + using field_type = typename V::field_type; + using GPUDILU = typename gpuistl::GpuDILU, gpuistl::GpuVector>; + return std::make_shared>(std::make_shared(op.getmat(), split_matrix, tune_gpu_kernels, mixed_precision_scheme)); + }); + + F::addCreator("GPUDILUFloat", [](const O& op, [[maybe_unused]] const P& prm, const std::function&, std::size_t) { + const bool split_matrix = prm.get("split_matrix", true); + const bool tune_gpu_kernels = prm.get("tune_gpu_kernels", true); + const int mixed_precision_scheme = prm.get("mixed_precision_scheme", 0); + + using block_type = typename V::block_type; + using VTo = Dune::BlockVector>; + using matrix_type_to = typename Dune::BCRSMatrix>; + using GpuDILU = typename gpuistl::GpuDILU, gpuistl::GpuVector>; + using Adapter = typename gpuistl::PreconditionerAdapter; + using Converter = typename gpuistl::PreconditionerConvertFieldTypeAdapter; + auto converted = std::make_shared(op.getmat()); + auto adapted = std::make_shared(std::make_shared(converted->getConvertedMatrix(), split_matrix, tune_gpu_kernels, mixed_precision_scheme)); + converted->setUnderlyingPreconditioner(adapted); + return converted; + }); +#endif + } +}; + + +} // namespace Opm + + +#endif // OPM_STANDARDPRECONDITIONERS_SERIAL_HPP \ No newline at end of file diff --git a/opm/simulators/linalg/is_gpu_operator.hpp b/opm/simulators/linalg/is_gpu_operator.hpp new file mode 100644 index 0000000000..a0a416ba5f --- /dev/null +++ b/opm/simulators/linalg/is_gpu_operator.hpp @@ -0,0 +1,55 @@ +/* + Copyright 2025 Equinor ASA + + This file is part of the Open Porous Media project (OPM). + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_IS_GPU_OPERATOR_HEADER +#define OPM_IS_GPU_OPERATOR_HEADER + +#include + +#include +#if HAVE_CUDA +#include +#endif +namespace Opm +{ + +/** + * \brief Check if a given operator is a GPU operator. + * + * This is used to check if the operator is a GPU operator, which is used + * in say the preconditioner factory to specialize the StandardPreconditioners class + * + * \tparam T The type of the operator to check. + */ +template +struct is_gpu_operator { +#if HAVE_CUDA + // TODO: This can be done more thoroughly by checking if range and matrix also is a gpu operator, but this + // works for now. + static constexpr bool value + = std::is_same_v>; +#else + // If CUDA is not enabled, we assume that the operator is not a GPU operator. + static constexpr bool value = false; +#endif +}; + +template +static constexpr bool is_gpu_operator_v = is_gpu_operator::value; + +} // namespace Opm + +#endif diff --git a/tests/gpuistl/test_preconditioner_factory_gpu.cpp b/tests/gpuistl/test_preconditioner_factory_gpu.cpp new file mode 100644 index 0000000000..f815857c5f --- /dev/null +++ b/tests/gpuistl/test_preconditioner_factory_gpu.cpp @@ -0,0 +1,80 @@ +/* + Copyright 2025 Equinor ASA + + This file is part of the Open Porous Media project (OPM). + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include +#include +#include +#include + +#define BOOST_TEST_MODULE TestPreconditionerFactoryGPU + +#include +#include +#include + +#include +#include + +using Scalar = double; +constexpr static int Dim = 3; +using MatrixTypeCPU = Dune::BCRSMatrix>; +// using VectorType = Dune::BlockVector>; +using CommSeq = Dune::Amg::SequentialInformation; +// using GpuOperator = Dune::MatrixAdapter; + +using GpuMatrixType = Opm::gpuistl::GpuSparseMatrix; +using GpuVectorType = Opm::gpuistl::GpuVector; +using GpuOperatorType = Dune::MatrixAdapter; +using FactoryTypeGpu = Opm::PreconditionerFactory; + +BOOST_AUTO_TEST_CASE(TestMatrixAdapter) +{ + const std::size_t N = 10; + const std::size_t nonZeroes = N * 3 - 2; + MatrixTypeCPU matrix(N, N, nonZeroes, MatrixTypeCPU::row_wise); + for (auto row = matrix.createbegin(); row != matrix.createend(); ++row) { + // Add nonzeros for left neighbour, diagonal and right neighbour + if (row.index() > 0) { + row.insert(row.index() - 1); + } + row.insert(row.index()); + if (row.index() < matrix.N() - 1) { + row.insert(row.index() + 1); + } + } + + for (std::size_t i = 0; i < N; ++i) { + matrix[i][i] = -2; + if (i < int(N) - 1) { + matrix[i][i + 1] = 1; + } + + if (i > 0) { + matrix[i][i - 1] = 1; + } + } + GpuMatrixType gpuMatrix = GpuMatrixType::fromMatrix(matrix); + Dune::MatrixAdapter gpuOp(gpuMatrix); + + Opm::PropertyTree prm; + prm.put(std::string("preconditioner"), std::string("DILU")); + auto weightsCalculator = []() { return GpuVectorType(10); }; + + // TODO: Once the GPU preconditioners are implemented, we can remove the check throw + // and actually test the preconditioner creation. + BOOST_CHECK_THROW(FactoryTypeGpu::create(gpuOp, prm, weightsCalculator), std::invalid_argument); +}