diff --git a/opm/simulators/linalg/PreconditionerFactory.hpp b/opm/simulators/linalg/PreconditionerFactory.hpp index e8015cbc8..338f5f873 100644 --- a/opm/simulators/linalg/PreconditionerFactory.hpp +++ b/opm/simulators/linalg/PreconditionerFactory.hpp @@ -197,6 +197,46 @@ private: } } + /// 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 preconditiner. + template + static size_t interiorIfGhostLast(const CommArg& comm) + { + size_t interior_count = 0; + size_t highest_interior_index = 0; + const auto& is = comm.indexSet(); + for (const auto& ind : is) { + if (CommArg::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(); + } + } + + static PrecPtr + createParILU(const Operator& op, const boost::property_tree::ptree& prm, const Comm& comm, const int ilulevel) + { + 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 size_t num_interior = interiorIfGhostLast(comm); + return std::make_shared>( + op.getmat(), comm, w, Opm::MILU_VARIANT::ILU, num_interior, redblack, reorder_spheres); + } else { + return std::make_shared>( + op.getmat(), comm, ilulevel, w, Opm::MILU_VARIANT::ILU, redblack, reorder_spheres); + } + } + // Add a useful default set of preconditioners to the factory. // This is the default template, used for parallel preconditioners. // (Serial specialization below). @@ -210,22 +250,13 @@ private: using P = boost::property_tree::ptree; using C = Comm; doAddCreator("ILU0", [](const O& op, const P& prm, const std::function&, const C& comm) { - const double w = prm.get("relaxation", 1.0); - return std::make_shared>( - op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU); + return createParILU(op, prm, comm, 0); }); doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function&, const C& comm) { - const double w = prm.get("relaxation", 1.0); - const int n = prm.get("ilulevel", 0); - // Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner. - return std::make_shared>( - op.getmat(), comm, n, w, Opm::MILU_VARIANT::ILU); + return createParILU(op, prm, comm, prm.get("ilulevel", 0)); }); doAddCreator("ILUn", [](const O& op, const P& prm, const std::function&, const C& comm) { - const int n = prm.get("ilulevel", 0); - const double w = prm.get("relaxation", 1.0); - return std::make_shared>( - op.getmat(), comm, n, w, Opm::MILU_VARIANT::ILU); + return createParILU(op, prm, comm, prm.get("ilulevel", 0)); }); doAddCreator("Jac", [](const O& op, const P& prm, const std::function&, const C& comm) { diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 7ed6c19f4..1c4b0db2f 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -26,7 +26,6 @@ #include #include -#include #include #include diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index c9c92a303..095c89a07 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -22,6 +22,7 @@ #include #include #include +#include namespace Opm {