19#ifndef OPM_SOLVERADAPTER_HPP
20#define OPM_SOLVERADAPTER_HPP
22#include <dune/istl/bcrsmatrix.hh>
23#include <dune/istl/operators.hh>
24#include <dune/istl/paamg/pinfo.hh>
25#include <dune/istl/schwarz.hh>
26#include <dune/istl/solver.hh>
28#include <opm/common/ErrorMacros.hpp>
55template <
class Operator,
template <
class>
class UnderlyingSolver,
class X>
59 using typename Dune::IterativeSolver<X, X>::domain_type;
60 using typename Dune::IterativeSolver<X, X>::range_type;
61 using typename Dune::IterativeSolver<X, X>::field_type;
62 using typename Dune::IterativeSolver<X, X>::real_type;
63 using typename Dune::IterativeSolver<X, X>::scalar_real_type;
64 static constexpr auto block_size = domain_type::block_type::dimension;
65 using XGPU = Opm::gpuistl::GpuVector<real_type>;
80 Dune::ScalarProduct<X>& sp,
81 std::shared_ptr<Dune::Preconditioner<X, X>> prec,
82 scalar_real_type reduction,
86 :
Dune::IterativeSolver<X, X>(op, sp, *prec, reduction, maxit, verbose)
87 , m_opOnCPUWithMatrix(op)
89 , m_underlyingSolver(constructSolver(prec, reduction, maxit, verbose, comm))
93 "Currently we only support operators of type MatrixAdapter in the CUDA/HIP solver. "
94 "Use --matrix-add-well-contributions=true. "
95 "Using WellModelMatrixAdapter with SolverAdapter is not well-defined so it will not work well, or at all.");
104 if (!m_inputBuffer) {
105 m_inputBuffer.reset(
new XGPU(b.dim()));
106 m_outputBuffer.reset(
new XGPU(x.dim()));
109 m_inputBuffer->copyFromHost(b);
111 m_outputBuffer->copyFromHost(x);
113 m_underlyingSolver.apply(*m_outputBuffer, *m_inputBuffer, reduction, res);
116 m_inputBuffer->copyToHost(b);
117 m_outputBuffer->copyToHost(x);
124 if (!m_inputBuffer) {
125 m_inputBuffer.reset(
new XGPU(b.dim()));
126 m_outputBuffer.reset(
new XGPU(x.dim()));
129 m_inputBuffer->copyFromHost(b);
131 m_outputBuffer->copyFromHost(x);
133 m_underlyingSolver.apply(*m_outputBuffer, *m_inputBuffer, res);
136 m_inputBuffer->copyToHost(b);
137 m_outputBuffer->copyToHost(x);
141 Operator& m_opOnCPUWithMatrix;
144 UnderlyingSolver<XGPU> m_underlyingSolver;
150 template <
class Comm>
151 UnderlyingSolver<XGPU> constructSolver(std::shared_ptr<Dune::Preconditioner<X, X>> prec,
152 scalar_real_type reduction,
155 const Comm& communication)
161 auto precAsHolder = std::dynamic_pointer_cast<PreconditionerHolder<X, X>>(prec);
163 OPM_THROW(std::invalid_argument,
164 "The preconditioner needs to be a CUDA preconditioner (eg. GPUDILU) wrapped in a "
165 "Opm::gpuistl::PreconditionerAdapter wrapped in a "
166 "Opm::gpuistl::GpuBlockPreconditioner. If you are unsure what this means, set "
167 "preconditioner to 'GPUDILU'");
170 auto preconditionerAdapter = precAsHolder->getUnderlyingPreconditioner();
171 auto preconditionerAdapterAsHolder
172 = std::dynamic_pointer_cast<PreconditionerHolder<XGPU, XGPU>>(preconditionerAdapter);
173 if (!preconditionerAdapterAsHolder) {
174 OPM_THROW(std::invalid_argument,
175 "The preconditioner needs to be a CUDA preconditioner (eg. GPUDILU) wrapped in a "
176 "Opm::gpuistl::PreconditionerAdapter wrapped in a "
177 "Opm::gpuistl::GpuBlockPreconditioner. If you are unsure what this means, set "
178 "preconditioner to 'GPUDILU'");
181 auto preconditionerReallyOnGPU = preconditionerAdapterAsHolder->getUnderlyingPreconditioner();
184 bool mpiSupportsCudaAwareAtCompileTime =
false;
185 bool mpiSupportsCudaAwareAtRunTime =
false;
187#if defined(MPIX_CUDA_AWARE_SUPPORT) && MPIX_CUDA_AWARE_SUPPORT
188 mpiSupportsCudaAwareAtCompileTime =
true;
191#if defined(MPIX_CUDA_AWARE_SUPPORT)
192 if (1 == MPIX_Query_cuda_support()) {
193 mpiSupportsCudaAwareAtRunTime =
true;
199 std::shared_ptr<Opm::gpuistl::GPUSender<real_type, Comm>> gpuComm;
200 if (mpiSupportsCudaAwareAtCompileTime && mpiSupportsCudaAwareAtRunTime) {
201 gpuComm = std::make_shared<
205 gpuComm = std::make_shared<
210 using CudaCommunication = GpuOwnerOverlapCopy<real_type, block_size, Comm>;
211 using SchwarzOperator
212 = Dune::OverlappingSchwarzOperator<GpuSparseMatrix<real_type>,
XGPU,
XGPU, CudaCommunication>;
213 auto cudaCommunication = std::make_shared<CudaCommunication>(gpuComm);
215 auto mpiPreconditioner = std::make_shared<GpuBlockPreconditioner<XGPU, XGPU, CudaCommunication>>(
216 preconditionerReallyOnGPU, cudaCommunication);
218 auto scalarProduct = std::make_shared<Dune::ParallelScalarProduct<XGPU, CudaCommunication>>(
219 cudaCommunication, m_opOnCPUWithMatrix.category());
226 OPM_ERROR_IF(cudaCommunication.use_count() < 2,
"Internal error. Shared pointer not owned properly.");
227 auto overlappingCudaOperator = std::make_shared<SchwarzOperator>(m_matrix, *cudaCommunication);
229 return UnderlyingSolver<XGPU>(
230 overlappingCudaOperator, scalarProduct, mpiPreconditioner, reduction, maxit, verbose);
235 UnderlyingSolver<XGPU> constructSolver(std::shared_ptr<Dune::Preconditioner<X, X>> prec,
236 scalar_real_type reduction,
239 [[maybe_unused]]
const Dune::Amg::SequentialInformation& communication)
242 return constructSolver(prec, reduction, maxit, verbose);
246 UnderlyingSolver<XGPU> constructSolver(std::shared_ptr<Dune::Preconditioner<X, X>> prec,
247 scalar_real_type reduction,
255 auto precAsHolder = std::dynamic_pointer_cast<PreconditionerHolder<XGPU, XGPU>>(prec);
257 OPM_THROW(std::invalid_argument,
258 "The preconditioner needs to be a CUDA preconditioner wrapped in a "
259 "Opm::gpuistl::PreconditionerHolder (eg. GPUDILU).");
261 auto preconditionerOnGPU = precAsHolder->getUnderlyingPreconditioner();
263 auto matrixOperator = std::make_shared<Dune::MatrixAdapter<GpuSparseMatrix<real_type>,
XGPU,
XGPU>>(m_matrix);
264 auto scalarProduct = std::make_shared<Dune::SeqScalarProduct<XGPU>>();
265 return UnderlyingSolver<XGPU>(matrixOperator, scalarProduct, preconditionerOnGPU, reduction, maxit, verbose);
268 std::unique_ptr<XGPU> m_inputBuffer;
269 std::unique_ptr<XGPU> m_outputBuffer;
Dune::OwnerOverlapCopyCommunication< int, int > Comm
Definition: FlexibleSolver_impl.hpp:304
Derived class of GPUSender that handles MPI made with CUDA aware MPI The copOwnerToAll function uses ...
Definition: GpuOwnerOverlapCopy.hpp:177
Derived class of GPUSender that handles MPI calls that should NOT use GPU direct communicatoin The im...
Definition: GpuOwnerOverlapCopy.hpp:121
The GpuSparseMatrix class simple wrapper class for a CuSparse matrix.
Definition: GpuSparseMatrix.hpp:48
void updateNonzeroValues(const MatrixType &matrix, bool copyNonZeroElementsDirectly=false)
updateNonzeroValues updates the non-zero values by using the non-zero values of the supplied matrix
Wraps a CUDA solver to work with CPU data.
Definition: SolverAdapter.hpp:57
virtual void apply(X &x, X &b, Dune::InverseOperatorResult &res) override
Definition: SolverAdapter.hpp:119
virtual void apply(X &x, X &b, double reduction, Dune::InverseOperatorResult &res) override
Definition: SolverAdapter.hpp:98
Opm::gpuistl::GpuVector< real_type > XGPU
Definition: SolverAdapter.hpp:65
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:79
static constexpr auto block_size
Definition: SolverAdapter.hpp:64
The is_a_well_operator class tries to guess if the operator is a well type operator.
Definition: has_function.hpp:114
Definition: fvbaseprimaryvariables.hh:141
Definition: autotuner.hpp:30
Dune::InverseOperatorResult InverseOperatorResult
Definition: GpuBridge.hpp:32