Loading...
Searching...
No Matches
FlexibleSolver_impl.hpp
/*
Copyright 2019, 2020 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2020 Equinor.
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
#define OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/TimingMacros.hpp>
#include <opm/simulators/linalg/matrixblock.hh>
#include <opm/simulators/linalg/ilufirstelement.hh>
#include <opm/simulators/linalg/FlexibleSolver.hpp>
#include <opm/simulators/linalg/PreconditionerFactory.hpp>
#include <opm/simulators/linalg/PropertyTree.hpp>
#include <opm/simulators/linalg/WellOperators.hpp>
#include <opm/simulators/linalg/PreconditionerFactoryGPUIncludeWrapper.hpp>
#include <opm/simulators/linalg/is_gpu_operator.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/umfpack.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <type_traits>
#if HAVE_CUDA
#if USE_HIP
#include <opm/simulators/linalg/gpuistl_hip/SolverAdapter.hpp>
#else
#include <opm/simulators/linalg/gpuistl/SolverAdapter.hpp>
#endif
#endif
namespace Dune
{
template <class Operator>
FlexibleSolver<Operator>::
FlexibleSolver(Operator& op,
const Opm::PropertyTree& prm,
const std::function<VectorType()>& weightsCalculator,
std::size_t pressureIndex)
{
init(op, Dune::Amg::SequentialInformation(), prm, weightsCalculator,
pressureIndex);
}
template <class Operator>
template <class Comm>
FlexibleSolver<Operator>::
FlexibleSolver(Operator& op,
const Comm& comm,
const Opm::PropertyTree& prm,
const std::function<VectorType()>& weightsCalculator,
std::size_t pressureIndex)
{
init(op, comm, prm, weightsCalculator, pressureIndex);
}
template <class Operator>
void
FlexibleSolver<Operator>::
apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res)
{
if (direct_solver_) {
recreateDirectSolver();
}
linsolver_->apply(x, rhs, res);
}
template <class Operator>
void
FlexibleSolver<Operator>::
apply(VectorType& x, VectorType& rhs, double reduction, Dune::InverseOperatorResult& res)
{
if (direct_solver_) {
recreateDirectSolver();
}
linsolver_->apply(x, rhs, reduction, res);
}
template <class Operator>
auto
FlexibleSolver<Operator>::
preconditioner() -> AbstractPrecondType&
{
return *preconditioner_;
}
template <class Operator>
Dune::SolverCategory::Category
FlexibleSolver<Operator>::
category() const
{
return linearoperator_for_solver_->category();
}
// Machinery for making sequential or parallel operators/preconditioners/scalar products.
template <class Operator>
template <class Comm>
void
FlexibleSolver<Operator>::
initOpPrecSp(Operator& op,
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator,
const Comm& comm,
std::size_t pressureIndex)
{
// Parallel case.
linearoperator_for_solver_ = &op;
auto child = prm.get_child_optional("preconditioner");
preconditioner_ = Opm::PreconditionerFactory<Operator, Comm>::create(op,
child ? *child : Opm::PropertyTree(),
weightsCalculator,
comm,
pressureIndex);
scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
}
template <class Operator>
void
FlexibleSolver<Operator>::
initOpPrecSp(Operator& op,
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator,
const Dune::Amg::SequentialInformation&,
std::size_t pressureIndex)
{
// Sequential case.
linearoperator_for_solver_ = &op;
auto child = prm.get_child_optional("preconditioner");
preconditioner_ = Opm::PreconditionerFactory<Operator,Dune::Amg::SequentialInformation>::create(op,
child ? *child : Opm::PropertyTree(),
weightsCalculator,
pressureIndex);
scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
}
template <class Operator>
template <class Comm>
void
FlexibleSolver<Operator>::
initSolver(const Opm::PropertyTree& prm, const Comm& comm)
{
const bool is_iorank = comm.communicator().rank() == 0;
const double tol = prm.get<double>("tol", 1e-2);
const int maxiter = prm.get<int>("maxiter", 200);
const int verbosity = is_iorank ? prm.get<int>("verbosity", 0) : 0;
const std::string solver_type = prm.get<std::string>("solver", "bicgstab");
// make sure it is nullptr at the start (used for error checking in the end).
// while the linSolver_ is initalized as a nullptr, we want to make sure it is reset here,
// simply because we will check if it is at the end of this function and need to keep this invariant
// (that it is nullptr at the start of this function).
linsolver_.reset();
if (solver_type == "bicgstab") {
linsolver_ = std::make_shared<Dune::BiCGSTABSolver<VectorType>>(*linearoperator_for_solver_,
*scalarproduct_,
*preconditioner_,
tol, // desired residual reduction factor
maxiter, // maximum number of iterations
verbosity);
} else if (solver_type == "loopsolver") {
linsolver_ = std::make_shared<Dune::LoopSolver<VectorType>>(*linearoperator_for_solver_,
*scalarproduct_,
*preconditioner_,
tol, // desired residual reduction factor
maxiter, // maximum number of iterations
verbosity);
} else if (solver_type == "gmres") {
int restart = prm.get<int>("restart", 15);
linsolver_ = std::make_shared<Dune::RestartedGMResSolver<VectorType>>(*linearoperator_for_solver_,
*scalarproduct_,
*preconditioner_,
tol,// desired residual reduction factor
restart,
maxiter, // maximum number of iterations
verbosity);
} else {
if constexpr (!Opm::is_gpu_operator_v<Operator>) {
if (solver_type == "flexgmres") {
int restart = prm.get<int>("restart", 15);
linsolver_ = std::make_shared<Dune::RestartedFlexibleGMResSolver<VectorType>>(*linearoperator_for_solver_,
*scalarproduct_,
*preconditioner_,
tol,// desired residual reduction factor
restart,
maxiter, // maximum number of iterations
verbosity);
#if HAVE_SUITESPARSE_UMFPACK
} else if (solver_type == "umfpack") {
if constexpr (std::is_same_v<typename VectorType::field_type,float>) {
OPM_THROW(std::invalid_argument, "UMFPack cannot be used with floats");
} else {
using MatrixType = std::remove_const_t<std::remove_reference_t<decltype(linearoperator_for_solver_->getmat())>>;
linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), verbosity, false);
direct_solver_ = true;
}
#endif
#if HAVE_CUDA
} else if (solver_type == "gpubicgstab") {
linsolver_.reset(new Opm::gpuistl::SolverAdapter<Operator, Dune::BiCGSTABSolver, VectorType>(
*linearoperator_for_solver_,
*scalarproduct_,
preconditioner_,
tol, // desired residual reduction factor
maxiter, // maximum number of iterations
verbosity,
comm));
#endif
}
}
}
if (!linsolver_) {
OPM_THROW(std::invalid_argument,
"Properties: Solver " + solver_type + " not known.");
}
}
// For now, the only direct solver we support is UMFPACK from SuiteSparse.
// When the matrix is updated (keeping sparsity pattern) it is possible to
// exploit separation of symbolic and numeric factorization, but we do not
// do so at this point. For complete generality, the solver abstract class
// Dune::InverseOperator<> should be extended with an update() function.
template <class Operator>
void
FlexibleSolver<Operator>::
recreateDirectSolver()
{
#if HAVE_SUITESPARSE_UMFPACK
if constexpr (!Opm::is_gpu_operator_v<Operator>) {
if constexpr (std::is_same_v<typename VectorType::field_type, float>) {
OPM_THROW(std::invalid_argument, "UMFPack cannot be used with floats");
} else {
using MatrixType = std::remove_const_t<std::remove_reference_t<decltype(linearoperator_for_solver_->getmat())>>;
linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), 0, false);
}
}
#else
OPM_THROW(std::logic_error, "Direct solver specified, but the FlexibleSolver class was not compiled with SuiteSparse support.");
#endif
}
// Main initialization routine.
// Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
template <class Operator>
template <class Comm>
void
FlexibleSolver<Operator>::
init(Operator& op,
const Comm& comm,
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator,
std::size_t pressureIndex)
{
initOpPrecSp(op, prm, weightsCalculator, comm, pressureIndex);
initSolver(prm, comm);
}
} // namespace Dune
// Macros to simplify explicit instantiation of FlexibleSolver for various block sizes.
// Vectors and matrices.
template<class Scalar, int N>
using BV = Dune::BlockVector<Dune::FieldVector<Scalar, N>>;
template<class Scalar, int N>
using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, N, N>>;
// Sequential operators.
template<class Scalar, int N>
using SeqOpM = Dune::MatrixAdapter<OBM<Scalar,N>, BV<Scalar,N>, BV<Scalar,N>>;
template<class Scalar, int N>
using SeqOpW = Opm::WellModelMatrixAdapter<OBM<Scalar,N>, BV<Scalar,N>, BV<Scalar,N>>;
#if HAVE_MPI
// Parallel communicator and operators.
using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
template<class Scalar, int N>
using ParOpM = Opm::GhostLastMatrixAdapter<OBM<Scalar,N>, BV<Scalar,N>, BV<Scalar,N>, Comm>;
template<class Scalar, int N>
using ParOpW = Opm::WellModelGhostLastMatrixAdapter<OBM<Scalar,N>, BV<Scalar,N>, BV<Scalar,N>, true>;
template<class Scalar, int N>
using ParOpD = Dune::OverlappingSchwarzOperator<OBM<Scalar,N>, BV<Scalar,N>, BV<Scalar,N>, Comm>;
// Note: we must instantiate the constructor that is a template.
// This is only needed in the parallel case, since otherwise the Comm type is
// not a template argument but always SequentialInformation.
#define INSTANTIATE_FLEXIBLESOLVER_OP(...) \
template class Dune::FlexibleSolver<__VA_ARGS__>; \
template Dune::FlexibleSolver<__VA_ARGS__>:: \
FlexibleSolver(__VA_ARGS__& op, \
const Comm& comm, \
const Opm::PropertyTree& prm, \
const std::function<typename __VA_ARGS__::domain_type()>& weightsCalculator, \
std::size_t pressureIndex);
#define INSTANTIATE_FLEXIBLESOLVER(T,N) \
INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<T,N>); \
INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<T,N>); \
INSTANTIATE_FLEXIBLESOLVER_OP(ParOpM<T,N>); \
INSTANTIATE_FLEXIBLESOLVER_OP(ParOpW<T,N>); \
INSTANTIATE_FLEXIBLESOLVER_OP(ParOpD<T,N>);
#else // HAVE_MPI
#define INSTANTIATE_FLEXIBLESOLVER_OP(...) \
template class Dune::FlexibleSolver<__VA_ARGS__>;
#define INSTANTIATE_FLEXIBLESOLVER(T,N) \
INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<T,N>); \
INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<T,N>);
#endif // HAVE_MPI
#endif // OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
Dune::FlexibleSolver
A solver class that encapsulates all needed objects for a linear solver (operator,...
Definition FlexibleSolver.hpp:45
Dune::FlexibleSolver::FlexibleSolver
FlexibleSolver(Operator &op, const Opm::PropertyTree &prm, const std::function< VectorType()> &weightsCalculator, std::size_t pressureIndex)
Create a sequential solver.
Definition FlexibleSolver_impl.hpp:56
Dune::FlexibleSolver::AbstractPrecondType
Dune::PreconditionerWithUpdate< VectorType, VectorType > AbstractPrecondType
Base class type of the contained preconditioner.
Definition FlexibleSolver.hpp:50
Dune::FlexibleSolver::preconditioner
AbstractPrecondType & preconditioner()
Access the contained preconditioner.
Definition FlexibleSolver_impl.hpp:104
Opm::GhostLastMatrixAdapter
Dune linear operator that assumes ghost rows are ordered after interior rows.
Definition WellOperators.hpp:402
Opm::PreconditionerFactory::create
static PrecPtr create(const Operator &op, const PropertyTree &prm, const std::function< Vector()> &weightsCalculator={}, std::size_t pressureIndex=std::numeric_limits< std::size_t >::max())
Create a new serial preconditioner and return a pointer to it.
Definition PreconditionerFactory_impl.hpp:154
Opm::PropertyTree
Hierarchical collection of key/value pairs.
Definition PropertyTree.hpp:39
Opm::PropertyTree::get_child_optional
std::optional< PropertyTree > get_child_optional(const std::string &key) const
Retrieve copy of sub tree rooted at node.
Definition PropertyTree.cpp:90
Opm::PropertyTree::get
T get(const std::string &key) const
Retrieve property value given hierarchical property key.
Definition PropertyTree.cpp:59
Opm::WellModelGhostLastMatrixAdapter
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:299
Opm::WellModelMatrixAdapter
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:225