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>
50template <
class Operator,
template <
class>
class UnderlyingSolver,
class X>
54 using typename Dune::IterativeSolver<X, X>::domain_type;
55 using typename Dune::IterativeSolver<X, X>::range_type;
56 using typename Dune::IterativeSolver<X, X>::field_type;
57 using typename Dune::IterativeSolver<X, X>::real_type;
58 using typename Dune::IterativeSolver<X, X>::scalar_real_type;
59 static constexpr auto block_size = domain_type::block_type::dimension;
76 Dune::ScalarProduct<X>& sp,
77 std::shared_ptr<Dune::Preconditioner<X, X>> prec,
78 scalar_real_type reduction,
82 :
Dune::IterativeSolver<X, X>(op, sp, *prec, reduction, maxit, verbose)
83 , m_opOnCPUWithMatrix(op)
85 , m_underlyingSolver(constructSolver(prec, reduction, maxit, verbose, comm))
89 "Currently we only support operators of type MatrixAdapter in the CUDA/HIP solver. "
90 "Use --matrix-add-well-contributions=true. "
91 "Using WellModelMatrixAdapter with SolverAdapter is not well-defined so it will not work well, or at all.");
100 if (!m_inputBuffer) {
101 m_inputBuffer.reset(
new XGPU(b.dim()));
102 m_outputBuffer.reset(
new XGPU(x.dim()));
105 m_inputBuffer->copyFromHost(b);
107 m_outputBuffer->copyFromHost(x);
109 m_underlyingSolver.apply(*m_outputBuffer, *m_inputBuffer, reduction, res);
112 m_inputBuffer->copyToHost(b);
113 m_outputBuffer->copyToHost(x);
120 if (!m_inputBuffer) {
121 m_inputBuffer.reset(
new XGPU(b.dim()));
122 m_outputBuffer.reset(
new XGPU(x.dim()));
125 m_inputBuffer->copyFromHost(b);
127 m_outputBuffer->copyFromHost(x);
129 m_underlyingSolver.apply(*m_outputBuffer, *m_inputBuffer, res);
132 m_inputBuffer->copyToHost(b);
133 m_outputBuffer->copyToHost(x);
137 Operator& m_opOnCPUWithMatrix;
140 UnderlyingSolver<XGPU> m_underlyingSolver;
146 template <
class Comm>
147 UnderlyingSolver<XGPU> constructSolver(std::shared_ptr<Dune::Preconditioner<X, X>> prec,
148 scalar_real_type reduction,
151 const Comm& communication)
157 auto precAsHolder = std::dynamic_pointer_cast<PreconditionerHolder<X, X>>(prec);
159 OPM_THROW(std::invalid_argument,
160 "The preconditioner needs to be a CUDA preconditioner (eg. GPUDILU) wrapped in a "
161 "Opm::gpuistl::PreconditionerAdapter wrapped in a "
162 "Opm::gpuistl::GpuBlockPreconditioner. If you are unsure what this means, set "
163 "preconditioner to 'gpudilu'");
166 auto preconditionerAdapter = precAsHolder->getUnderlyingPreconditioner();
167 auto preconditionerAdapterAsHolder
168 = std::dynamic_pointer_cast<PreconditionerHolder<XGPU, XGPU>>(preconditionerAdapter);
169 if (!preconditionerAdapterAsHolder) {
170 OPM_THROW(std::invalid_argument,
171 "The preconditioner needs to be a CUDA preconditioner (eg. GPUDILU) wrapped in a "
172 "Opm::gpuistl::PreconditionerAdapter wrapped in a "
173 "Opm::gpuistl::GpuBlockPreconditioner. If you are unsure what this means, set "
174 "preconditioner to 'gpudilu'");
177 auto preconditionerReallyOnGPU = preconditionerAdapterAsHolder->getUnderlyingPreconditioner();
180 using CudaCommunication = GpuOwnerOverlapCopy<real_type, Comm>;
181 using SchwarzOperator
182 = Dune::OverlappingSchwarzOperator<GpuSparseMatrix<real_type>,
XGPU,
XGPU, CudaCommunication>;
183 auto cudaCommunication = makeGpuOwnerOverlapCopy<real_type, block_size, Comm>(communication);
185 auto mpiPreconditioner = std::make_shared<GpuBlockPreconditioner<XGPU, XGPU, CudaCommunication>>(
186 preconditionerReallyOnGPU, cudaCommunication);
188 auto scalarProduct = std::make_shared<Dune::ParallelScalarProduct<XGPU, CudaCommunication>>(
189 cudaCommunication, m_opOnCPUWithMatrix.category());
196 OPM_ERROR_IF(cudaCommunication.use_count() < 2,
"Internal error. Shared pointer not owned properly.");
197 auto overlappingCudaOperator = std::make_shared<SchwarzOperator>(m_matrix, *cudaCommunication);
199 return UnderlyingSolver<XGPU>(
200 overlappingCudaOperator, scalarProduct, mpiPreconditioner, reduction, maxit, verbose);
205 UnderlyingSolver<XGPU> constructSolver(std::shared_ptr<Dune::Preconditioner<X, X>> prec,
206 scalar_real_type reduction,
209 [[maybe_unused]]
const Dune::Amg::SequentialInformation& communication)
212 return constructSolver(prec, reduction, maxit, verbose);
216 UnderlyingSolver<XGPU> constructSolver(std::shared_ptr<Dune::Preconditioner<X, X>> prec,
217 scalar_real_type reduction,
225 auto precAsHolder = std::dynamic_pointer_cast<PreconditionerHolder<XGPU, XGPU>>(prec);
227 OPM_THROW(std::invalid_argument,
228 "The preconditioner needs to be a CUDA preconditioner wrapped in a "
229 "Opm::gpuistl::PreconditionerHolder (eg. GPUDILU).");
231 auto preconditionerOnGPU = precAsHolder->getUnderlyingPreconditioner();
233 auto matrixOperator = std::make_shared<Dune::MatrixAdapter<GpuSparseMatrix<real_type>,
XGPU,
XGPU>>(m_matrix);
234 auto scalarProduct = std::make_shared<Dune::SeqScalarProduct<XGPU>>();
235 return UnderlyingSolver<XGPU>(matrixOperator, scalarProduct, preconditionerOnGPU, reduction, maxit, verbose);
238 std::unique_ptr<XGPU> m_inputBuffer;
239 std::unique_ptr<XGPU> m_outputBuffer;
Dune::OwnerOverlapCopyCommunication< int, int > Comm
Definition: FlexibleSolver_impl.hpp:304
The GpuSparseMatrix class simple wrapper class for a CuSparse matrix.
Definition: GpuSparseMatrix.hpp:60
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:52
virtual void apply(X &x, X &b, Dune::InverseOperatorResult &res) override
Definition: SolverAdapter.hpp:115
virtual void apply(X &x, X &b, double reduction, Dune::InverseOperatorResult &res) override
Definition: SolverAdapter.hpp:94
Opm::gpuistl::GpuVector< real_type > XGPU
Definition: SolverAdapter.hpp:60
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
static constexpr auto block_size
Definition: SolverAdapter.hpp:59
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: AmgxInterface.hpp:38
Dune::InverseOperatorResult InverseOperatorResult
Definition: GpuBridge.hpp:32