19#ifndef OPM_SOLVERADAPTER_HPP
20#define OPM_SOLVERADAPTER_HPP
24#include <dune/istl/bcrsmatrix.hh>
25#include <dune/istl/operators.hh>
26#include <dune/istl/schwarz.hh>
27#include <dune/istl/solver.hh>
28#include <opm/common/ErrorMacros.hpp>
49template <
class Operator,
template <
class>
class UnderlyingSolver,
class X>
53 using typename Dune::IterativeSolver<X, X>::domain_type;
54 using typename Dune::IterativeSolver<X, X>::range_type;
55 using typename Dune::IterativeSolver<X, X>::field_type;
56 using typename Dune::IterativeSolver<X, X>::real_type;
57 using typename Dune::IterativeSolver<X, X>::scalar_real_type;
58 static constexpr auto block_size = domain_type::block_type::dimension;
59 using XGPU = Opm::cuistl::CuVector<real_type>;
63 Dune::ScalarProduct<X>& sp,
64 std::shared_ptr<Dune::Preconditioner<X, X>> prec,
65 scalar_real_type reduction,
68 :
Dune::IterativeSolver<X, X>(op, sp, *prec, reduction, maxit, verbose)
69 , m_opOnCPUWithMatrix(op)
71 , m_underlyingSolver(constructSolver(prec, reduction, maxit, verbose))
82 m_inputBuffer.reset(
new XGPU(b.dim()));
83 m_outputBuffer.reset(
new XGPU(x.dim()));
86 m_inputBuffer->copyFromHost(b);
88 m_outputBuffer->copyFromHost(x);
90 m_underlyingSolver.apply(*m_outputBuffer, *m_inputBuffer, reduction, res);
93 m_inputBuffer->copyToHost(b);
94 m_outputBuffer->copyToHost(x);
101 if (!m_inputBuffer) {
102 m_inputBuffer.reset(
new XGPU(b.dim()));
103 m_outputBuffer.reset(
new XGPU(x.dim()));
106 m_inputBuffer->copyFromHost(b);
108 m_outputBuffer->copyFromHost(x);
110 m_underlyingSolver.apply(*m_outputBuffer, *m_inputBuffer, res);
113 m_inputBuffer->copyToHost(b);
114 m_outputBuffer->copyToHost(x);
118 Operator& m_opOnCPUWithMatrix;
121 UnderlyingSolver<XGPU> m_underlyingSolver;
125 UnderlyingSolver<XGPU> constructSolver(std::shared_ptr<Dune::Preconditioner<X, X>> prec,
126 scalar_real_type reduction,
133 "Currently we only support operators of type MatrixAdapter in the CUDA solver. "
134 "Use --matrix-add-well-contributions=true. "
135 "Using WellModelMatrixAdapter with SolverAdapter is not well-defined so it will not work well, or at all.");
147 auto precAsHolder = std::dynamic_pointer_cast<PreconditionerHolder<X, X>>(prec);
149 OPM_THROW(std::invalid_argument,
150 "The preconditioner needs to be a CUDA preconditioner (eg. CuILU0) wrapped in a "
151 "Opm::cuistl::PreconditionerAdapter wrapped in a "
152 "Opm::cuistl::CuBlockPreconditioner. If you are unsure what this means, set "
153 "preconditioner to 'CUILU0'");
156 auto preconditionerAdapter = precAsHolder->getUnderlyingPreconditioner();
157 auto preconditionerAdapterAsHolder
158 = std::dynamic_pointer_cast<PreconditionerHolder<XGPU, XGPU>>(preconditionerAdapter);
159 if (!preconditionerAdapterAsHolder) {
160 OPM_THROW(std::invalid_argument,
161 "The preconditioner needs to be a CUDA preconditioner (eg. CuILU0) wrapped in a "
162 "Opm::cuistl::PreconditionerAdapter wrapped in a "
163 "Opm::cuistl::CuBlockPreconditioner. If you are unsure what this means, set "
164 "preconditioner to 'CUILU0'");
167 auto preconditionerReallyOnGPU = preconditionerAdapterAsHolder->getUnderlyingPreconditioner();
168 const auto& communication = m_opOnCPUWithMatrix.getCommunication();
171 bool mpiSUpportsCudaAwareAtCompileTime =
false;
172 bool mpiSupportsCudaAwareAtRunTime =
false;
174#if defined(MPIX_CUDA_AWARE_SUPPORT) && MPIX_CUDA_AWARE_SUPPORT
175 mpiSupportsCudaAwareAtCompileTime =
true;
178#if defined(MPIX_CUDA_AWARE_SUPPORT)
179 if (1 == MPIX_Query_cuda_support()) {
180 mpiSupportsCudaAwareAtRunTime =
true;
186 std::shared_ptr<Opm::cuistl::GPUSender<real_type, typename Operator::communication_type>> gpuComm;
187 if (mpiSUpportsCudaAwareAtCompileTime && mpiSupportsCudaAwareAtRunTime){
188 gpuComm = std::make_shared<Opm::cuistl::GPUAwareMPISender<real_type, block_size, typename Operator::communication_type>>(communication);
191 gpuComm = std::make_shared<Opm::cuistl::GPUObliviousMPISender<real_type, block_size, typename Operator::communication_type>>(communication);
194 using CudaCommunication = CuOwnerOverlapCopy<real_type, block_size, typename Operator::communication_type>;
195 using SchwarzOperator
196 = Dune::OverlappingSchwarzOperator<CuSparseMatrix<real_type>,
XGPU,
XGPU, CudaCommunication>;
197 auto cudaCommunication = std::make_shared<CudaCommunication>(gpuComm);
199 auto mpiPreconditioner = std::make_shared<CuBlockPreconditioner<XGPU, XGPU, CudaCommunication>>(
200 preconditionerReallyOnGPU, cudaCommunication);
202 auto scalarProduct = std::make_shared<Dune::ParallelScalarProduct<XGPU, CudaCommunication>>(
203 cudaCommunication, m_opOnCPUWithMatrix.category());
211 OPM_ERROR_IF(cudaCommunication.use_count() < 2,
"Internal error. Shared pointer not owned properly.");
212 auto overlappingCudaOperator = std::make_shared<SchwarzOperator>(m_matrix, *cudaCommunication);
214 return UnderlyingSolver<XGPU>(
215 overlappingCudaOperator, scalarProduct, mpiPreconditioner, reduction, maxit, verbose);
221 auto precAsHolder = std::dynamic_pointer_cast<PreconditionerHolder<XGPU, XGPU>>(prec);
223 OPM_THROW(std::invalid_argument,
224 "The preconditioner needs to be a CUDA preconditioner wrapped in a "
225 "Opm::cuistl::PreconditionerHolder (eg. CuILU0).");
227 auto preconditionerOnGPU = precAsHolder->getUnderlyingPreconditioner();
230 = std::make_shared<Dune::MatrixAdapter<CuSparseMatrix<real_type>,
XGPU,
XGPU>>(m_matrix);
231 auto scalarProduct = std::make_shared<Dune::SeqScalarProduct<XGPU>>();
232 return UnderlyingSolver<XGPU>(
233 matrixOperator, scalarProduct, preconditionerOnGPU, reduction, maxit, verbose);
237 std::unique_ptr<XGPU> m_inputBuffer;
238 std::unique_ptr<XGPU> m_outputBuffer;
The CuSparseMatrix class simple wrapper class for a CuSparse matrix.
Definition: CuSparseMatrix.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:51
Opm::cuistl::CuVector< real_type > XGPU
Definition: SolverAdapter.hpp:59
virtual void apply(X &x, X &b, double reduction, Dune::InverseOperatorResult &res) override
Definition: SolverAdapter.hpp:75
SolverAdapter(Operator &op, Dune::ScalarProduct< X > &sp, std::shared_ptr< Dune::Preconditioner< X, X > > prec, scalar_real_type reduction, int maxit, int verbose)
Definition: SolverAdapter.hpp:62
static constexpr auto block_size
Definition: SolverAdapter.hpp:58
virtual void apply(X &x, X &b, Dune::InverseOperatorResult &res) override
Definition: SolverAdapter.hpp:96
The has_communication class checks if the type has the member function getCommunication.
Definition: has_function.hpp:96
The is_a_well_operator class tries to guess if the operator is a well type operator.
Definition: has_function.hpp:114
Definition: SupportsFaceTag.hpp:27
Definition: CuBlockPreconditioner.hpp:29
Dune::InverseOperatorResult InverseOperatorResult
Definition: BdaBridge.hpp:32