Loading...
Searching...
No Matches
SolverAdapter.hpp
/*
Copyright 2022-2023 SINTEF AS
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_SOLVERADAPTER_HPP
#define OPM_SOLVERADAPTER_HPP
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/solver.hh>
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/linalg/gpuistl/GpuBlockPreconditioner.hpp>
#include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
#include <opm/simulators/linalg/gpuistl/PreconditionerAdapter.hpp>
#include <opm/simulators/linalg/gpuistl/detail/has_function.hpp>
#if HAVE_MPI
#include <opm/simulators/linalg/gpuistl/GpuOwnerOverlapCopy.hpp>
#endif
#include <memory>
namespace Opm::gpuistl
{
template <class Operator, template <class> class UnderlyingSolver, class X>
class SolverAdapter : public Dune::IterativeSolver<X, X>
{
public:
using typename Dune::IterativeSolver<X, X>::domain_type;
using typename Dune::IterativeSolver<X, X>::range_type;
using typename Dune::IterativeSolver<X, X>::field_type;
using typename Dune::IterativeSolver<X, X>::real_type;
using typename Dune::IterativeSolver<X, X>::scalar_real_type;
static constexpr auto block_size = domain_type::block_type::dimension;
using XGPU = Opm::gpuistl::GpuVector<real_type>;
template<class Comm>
SolverAdapter(Operator& op,
Dune::ScalarProduct<X>& sp,
std::shared_ptr<Dune::Preconditioner<X, X>> prec,
scalar_real_type reduction,
int maxit,
int verbose,
const Comm& comm)
: Dune::IterativeSolver<X, X>(op, sp, *prec, reduction, maxit, verbose)
, m_opOnCPUWithMatrix(op)
, m_matrix(GpuSparseMatrixWrapper<real_type>::fromMatrix(op.getmat()))
, m_underlyingSolver(constructSolver(prec, reduction, maxit, verbose, comm))
{
OPM_ERROR_IF(
detail::is_a_well_operator<Operator>::value,
"Currently we only support operators of type MatrixAdapter in the CUDA/HIP solver. "
"Use --matrix-add-well-contributions=true. "
"Using WellModelMatrixAdapter with SolverAdapter is not well-defined so it will not work well, or at all.");
}
virtual void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override
{
// TODO: Can we do this without reimplementing the other function?
// TODO: [perf] Do we need to update the matrix every time? Probably yes
m_matrix.updateNonzeroValues(m_opOnCPUWithMatrix.getmat(), true);
if (!m_inputBuffer) {
m_inputBuffer.reset(new XGPU(b.dim()));
m_outputBuffer.reset(new XGPU(x.dim()));
}
m_inputBuffer->copyFromHost(b);
// TODO: [perf] do we need to copy x here?
m_outputBuffer->copyFromHost(x);
m_underlyingSolver.apply(*m_outputBuffer, *m_inputBuffer, reduction, res);
// TODO: [perf] do we need to copy b here?
m_inputBuffer->copyToHost(b);
m_outputBuffer->copyToHost(x);
}
virtual void apply(X& x, X& b, Dune::InverseOperatorResult& res) override
{
// TODO: [perf] Do we need to update the matrix every time? Probably yes
m_matrix.updateNonzeroValues(m_opOnCPUWithMatrix.getmat(), true);
if (!m_inputBuffer) {
m_inputBuffer.reset(new XGPU(b.dim()));
m_outputBuffer.reset(new XGPU(x.dim()));
}
m_inputBuffer->copyFromHost(b);
// TODO: [perf] do we need to copy x here?
m_outputBuffer->copyFromHost(x);
m_underlyingSolver.apply(*m_outputBuffer, *m_inputBuffer, res);
// TODO: [perf] do we need to copy b here?
m_inputBuffer->copyToHost(b);
m_outputBuffer->copyToHost(x);
}
private:
Operator& m_opOnCPUWithMatrix;
GpuSparseMatrixWrapper<real_type> m_matrix;
UnderlyingSolver<XGPU> m_underlyingSolver;
// TODO: Use a std::forward
// This is the MPI parallel case (general communication object)
#if HAVE_MPI
template <class Comm>
UnderlyingSolver<XGPU> constructSolver(std::shared_ptr<Dune::Preconditioner<X, X>> prec,
scalar_real_type reduction,
int maxit,
int verbose,
const Comm& communication)
{
// TODO: See the below TODO over the definition of precHolder in the other overload of constructSolver
// TODO: We are currently double wrapping preconditioners in the preconditioner factory to be extra
// compatible with CPU. Probably a cleaner way eventually would be to do more modifications to the
// flexible solver to accomodate the pure GPU better.
auto precAsHolder = std::dynamic_pointer_cast<PreconditionerHolder<X, X>>(prec);
if (!precAsHolder) {
OPM_THROW(std::invalid_argument,
"The preconditioner needs to be a CUDA preconditioner (eg. GPUDILU) wrapped in a "
"Opm::gpuistl::PreconditionerAdapter wrapped in a "
"Opm::gpuistl::GpuBlockPreconditioner. If you are unsure what this means, set "
"preconditioner to 'gpudilu'");
}
auto preconditionerAdapter = precAsHolder->getUnderlyingPreconditioner();
auto preconditionerAdapterAsHolder
= std::dynamic_pointer_cast<PreconditionerHolder<XGPU, XGPU>>(preconditionerAdapter);
if (!preconditionerAdapterAsHolder) {
OPM_THROW(std::invalid_argument,
"The preconditioner needs to be a CUDA preconditioner (eg. GPUDILU) wrapped in a "
"Opm::gpuistl::PreconditionerAdapter wrapped in a "
"Opm::gpuistl::GpuBlockPreconditioner. If you are unsure what this means, set "
"preconditioner to 'gpudilu'");
}
// We need to get the underlying preconditioner:
auto preconditionerReallyOnGPU = preconditionerAdapterAsHolder->getUnderlyingPreconditioner();
using CudaCommunication = GpuOwnerOverlapCopy<real_type, Comm>;
using SchwarzOperator
= Dune::OverlappingSchwarzOperator<GpuSparseMatrixWrapper<real_type>, XGPU, XGPU, CudaCommunication>;
auto cudaCommunication = makeGpuOwnerOverlapCopy<real_type, block_size, Comm>(communication);
auto mpiPreconditioner = std::make_shared<GpuBlockPreconditioner<XGPU, XGPU, CudaCommunication>>(
preconditionerReallyOnGPU, cudaCommunication);
auto scalarProduct = std::make_shared<Dune::ParallelScalarProduct<XGPU, CudaCommunication>>(
cudaCommunication, m_opOnCPUWithMatrix.category());
// NOTE: Ownership of cudaCommunication is handled by mpiPreconditioner. However, just to make sure we
// remember this, we add this check. So remember that we hold one count in this scope, and one in the
// GpuBlockPreconditioner instance. We accommodate for the fact that it could be passed around in
// GpuBlockPreconditioner, hence we do not test for != 2
OPM_ERROR_IF(cudaCommunication.use_count() < 2, "Internal error. Shared pointer not owned properly.");
auto overlappingCudaOperator = std::make_shared<SchwarzOperator>(m_matrix, *cudaCommunication);
return UnderlyingSolver<XGPU>(
overlappingCudaOperator, scalarProduct, mpiPreconditioner, reduction, maxit, verbose);
}
#endif
// This is the serial case (specific overload for Dune::Amg::SequentialInformation)
UnderlyingSolver<XGPU> constructSolver(std::shared_ptr<Dune::Preconditioner<X, X>> prec,
scalar_real_type reduction,
int maxit,
int verbose,
[[maybe_unused]] const Dune::Amg::SequentialInformation& communication)
{
// Dune::Amg::SequentialInformation is the serial case
return constructSolver(prec, reduction, maxit, verbose);
}
// TODO: Use a std::forward
UnderlyingSolver<XGPU> constructSolver(std::shared_ptr<Dune::Preconditioner<X, X>> prec,
scalar_real_type reduction,
int maxit,
int verbose)
{
// TODO: Fix the reliance on casting here. This is a code smell to a certain degree, and assumes
// a certain setup beforehand. The only reason we do it this way is to minimize edits to the
// flexible solver. We could design it differently, but keep this for the time being until
// we figure out how we want to GPU-ify the rest of the system.
auto precAsHolder = std::dynamic_pointer_cast<PreconditionerHolder<XGPU, XGPU>>(prec);
if (!precAsHolder) {
OPM_THROW(std::invalid_argument,
"The preconditioner needs to be a CUDA preconditioner wrapped in a "
"Opm::gpuistl::PreconditionerHolder (eg. GPUDILU).");
}
auto preconditionerOnGPU = precAsHolder->getUnderlyingPreconditioner();
auto matrixOperator = std::make_shared<Dune::MatrixAdapter<GpuSparseMatrixWrapper<real_type>, XGPU, XGPU>>(m_matrix);
auto scalarProduct = std::make_shared<Dune::SeqScalarProduct<XGPU>>();
return UnderlyingSolver<XGPU>(matrixOperator, scalarProduct, preconditionerOnGPU, reduction, maxit, verbose);
}
std::unique_ptr<XGPU> m_inputBuffer;
std::unique_ptr<XGPU> m_outputBuffer;
};
} // namespace Opm::gpuistl
#endif
Opm::gpuistl::GpuSparseMatrixWrapper
The GpuSparseMatrixWrapper Checks CUDA/HIP version and dispatches a version either using the old or t...
Definition GpuSparseMatrixWrapper.hpp:61
Opm::gpuistl::GpuSparseMatrixWrapper::updateNonzeroValues
void updateNonzeroValues(const MatrixType &matrix, bool copyNonZeroElementsDirectly=false)
updateNonzeroValues updates the non-zero values by using the non-zero values of the supplied matrix
Definition GpuSparseMatrixWrapper.hpp:374
Opm::gpuistl::GpuVector
Definition gpu_type_detection.hpp:30
Opm::gpuistl::SolverAdapter::SolverAdapter
SolverAdapter(Operator &op, Dune::ScalarProduct< X > &sp, std::shared_ptr< Dune::Preconditioner< X, X > > prec, scalar_real_type reduction, int maxit, int verbose, const Comm &comm)
constructor
Definition SolverAdapter.hpp:75