Loading...
Searching...
No Matches
ISTLSolverGPUISTL.hpp
/*
Copyright 2025 Equinor ASA
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_ISTLSOLVERGPUISTL_HEADER_INCLUDED
#define OPM_ISTLSOLVERGPUISTL_HEADER_INCLUDED
#include <dune/istl/operators.hh>
#include <memory>
#include <optional>
#include <opm/grid/utility/ElementChunks.hpp>
#include <opm/simulators/linalg/AbstractISTLSolver.hpp>
#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
#include <opm/simulators/linalg/ISTLSolver.hpp>
#if USE_HIP
#include <opm/simulators/linalg/gpuistl_hip/GpuSparseMatrixWrapper.hpp>
#include <opm/simulators/linalg/gpuistl_hip/GpuVector.hpp>
#include <opm/simulators/linalg/gpuistl_hip/PinnedMemoryHolder.hpp>
#else
#include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
#include <opm/simulators/linalg/gpuistl/PinnedMemoryHolder.hpp>
#endif
#include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
#include <opm/simulators/linalg/gpuistl/detail/FlexibleSolverWrapper.hpp>
#include <opm/simulators/linalg/printlinearsolverparameter.hpp>
namespace Opm::gpuistl
{
template <class TypeTag>
class ISTLSolverGPUISTL : public AbstractISTLSolver<TypeTag>
{
public:
using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
using Matrix = typename SparseMatrixAdapter::IstlMatrix;
using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using ElementChunksType = Opm::ElementChunks<GridView, Dune::Partitions::All>;
using real_type = typename Vector::field_type;
using GPUMatrix = Opm::gpuistl::GpuSparseMatrixWrapper<real_type>;
using GPUVector = Opm::gpuistl::GpuVector<real_type>;
using GPUVectorInt = Opm::gpuistl::GpuVector<int>;
constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
#if HAVE_MPI
using CommunicationType = Dune::OwnerOverlapCopyCommunication<int, int>;
#else
using CommunicationType = Dune::Communication<int>;
#endif
using SolverType = Opm::gpuistl::detail::FlexibleSolverWrapper<GPUMatrix, GPUVector, CommunicationType>;
ISTLSolverGPUISTL(const Simulator& simulator,
const FlowLinearSolverParameters& parameters,
bool forceSerial = false)
: m_parameters(parameters)
, m_simulator(simulator)
, m_forceSerial(forceSerial)
, m_element_chunks(simulator.vanguard().gridView(), Dune::Partitions::all, ThreadManager::maxThreads())
{
#if HAVE_MPI
m_comm = std::make_shared<CommunicationType>(simulator.vanguard().grid().comm());
// Extract parallel grid information to populate index sets
extractParallelGridInformationToISTL(simulator.vanguard().grid(), m_parallelInformation);
// Set up element mapper manually
ElementMapper elemMapper(simulator.vanguard().gridView(), Dune::mcmgElementLayout());
// Overlap rows are needed for making overlap rows invalid in parallel mode
Opm::detail::findOverlapAndInterior(simulator.vanguard().grid(), elemMapper, m_overlapRows, m_interiorRows);
if (isParallel()) {
const std::size_t size = simulator.vanguard().grid().leafGridView().size(0);
// Copy parallel information to communication object (index sets and remote indices)
Opm::detail::copyParValues(m_parallelInformation, size, *m_comm);
}
#else
m_comm = std::make_shared<CommunicationType>(simulator.gridView().comm());
#endif
m_parameters.init(simulator.vanguard().eclState().getSimulationConfig().useCPR());
m_propertyTree = setupPropertyTree(m_parameters,
Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
Parameters::IsSet<Parameters::LinearSolverReduction>());
if (!Parameters::Get<Parameters::MatrixAddWellContributions>()) {
OPM_THROW(std::logic_error, "Well operators are currently not supported for the GPU backend. "
"Use --matrix-add-well-contributions=true to add well contributions to the matrix instead.");
}
Opm::detail::printLinearSolverParameters(m_parameters, m_propertyTree, simulator.gridView().comm());
}
explicit ISTLSolverGPUISTL(const Simulator& simulator)
: ISTLSolverGPUISTL(simulator, FlowLinearSolverParameters(), false)
{
}
void eraseMatrix() override
{
// Nothing, this is the same as the ISTLSolver
}
void setActiveSolver(int num) override
{
if (num != 0) {
OPM_THROW(std::logic_error, "Only one solver available for the GPU backend.");
}
}
int numAvailableSolvers() const override
{
return 1;
}
void prepare(const SparseMatrixAdapter& M, Vector& b) override
{
prepare(M.istlMatrix(), b);
}
void prepare(const Matrix& M, Vector& b) override
{
try {
if (isParallel() && !m_overlapRows.empty()) {
Opm::detail::makeOverlapRowsInvalid(const_cast<Matrix&>(M), m_overlapRows);
}
updateMatrix(M);
updateRhs(b);
}
OPM_CATCH_AND_RETHROW_AS_CRITICAL_ERROR("This is likely due to a faulty linear solver JSON specification. "
"Check for errors related to missing nodes.");
}
void setResidual(Vector&) override
{
// Should be handled in prepare() instead.
}
void getResidual(Vector& b) const override
{
if (!m_rhs) {
OPM_THROW(std::runtime_error, "m_rhs not initialized, prepare(matrix, rhs); needs to be called");
}
m_rhs->copyToHost(b);
}
void setMatrix(const SparseMatrixAdapter&) override
{
// Should be handled in prepare() instead.
}
bool solve(Vector& x) override
{
// TODO: Write matrix to disk if needed
Dune::InverseOperatorResult result;
if (!m_matrix) {
OPM_THROW(std::runtime_error, "m_matrix not initialized, prepare(matrix, rhs); needs to be called");
}
if (!m_rhs) {
OPM_THROW(std::runtime_error, "m_rhs not initialized, prepare(matrix, rhs); needs to be called");
}
if (!m_gpuSolver) {
OPM_THROW(std::runtime_error,
"m_gpuFlexibleSolver not initialized, prepare(matrix, rhs); needs to be called");
}
if (!m_x) {
m_x = std::make_unique<GPUVector>(x);
m_pinnedXMemory = std::make_unique<PinnedMemoryHolder<real_type>>(
const_cast<real_type*>(&x[0][0]),
x.dim()
);
} else {
// copy from host to device using main stream and asynchronous transfer
m_x->copyFromHostAsync(x);
}
m_gpuSolver->apply(*m_x, *m_rhs, result);
m_x->copyToHost(x);
++m_solveCount;
m_lastSeenIterations = result.iterations;
return checkConvergence(result);
}
int iterations() const override
{
return m_lastSeenIterations;
}
const CommunicationType* comm() const override
{
return m_comm.get();
}
bool isParallel() const
{
#if HAVE_MPI
return !m_forceSerial && m_comm->communicator().size() > 1;
#else
return false;
#endif
}
int getSolveCount() const override
{
return m_solveCount;
}
private:
bool checkConvergence(const Dune::InverseOperatorResult& result) const
{
return AbstractISTLSolver<TypeTag>::checkConvergence(result, m_parameters);
}
// Weights to make approximate pressure equations.
std::function<GPUVector&()> getWeightsCalculator()
{
std::function<GPUVector&()> weightsCalculator;
using namespace std::string_literals;
auto preconditionerType = m_propertyTree.get("preconditioner.type"s, "cpr"s);
// Make the preconditioner type lowercase for internal canonical representation
std::transform(preconditionerType.begin(), preconditionerType.end(), preconditionerType.begin(), ::tolower);
if (preconditionerType == "cpr" || preconditionerType == "cprt"
|| preconditionerType == "cprw" || preconditionerType == "cprwt") {
const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
const auto weightsType = m_propertyTree.get("preconditioner.weight_type"s, "quasiimpes"s);
if (weightsType == "quasiimpes") {
m_weights.emplace(m_matrix->N() * m_matrix->blockSize());
// Pre-compute diagonal indices once when setting up the calculator
auto diagonalIndices = Amg::precomputeDiagonalIndices(*m_matrix);
m_diagonalIndices.emplace(diagonalIndices);
if (transpose) {
weightsCalculator = [this]() -> GPUVector& {
Amg::getQuasiImpesWeights<real_type, true>(
*m_matrix, pressureIndex, *m_weights, *m_diagonalIndices);
return *m_weights;
};
} else {
weightsCalculator = [this]() -> GPUVector& {
Amg::getQuasiImpesWeights<real_type, false>(
*m_matrix, pressureIndex, *m_weights, *m_diagonalIndices);
return *m_weights;
};
}
} else if (weightsType == "trueimpes") {
// Create CPU vector for the weights and initialize GPU vector
m_cpuWeights.resize(m_matrix->N());
m_pinnedWeightsMemory = std::make_unique<PinnedMemoryHolder<real_type>>(
const_cast<real_type*>(&m_cpuWeights[0][0]), m_cpuWeights.dim());
m_weights.emplace(m_cpuWeights);
const bool enableThreadParallel = m_parameters.cpr_weights_thread_parallel_;
// CPU implementation wrapped for GPU
weightsCalculator = [this, enableThreadParallel]() -> GPUVector& {
// Use the CPU implementation to calculate the weights
ElementContext elemCtx(m_simulator);
Amg::getTrueImpesWeights(pressureIndex,
m_cpuWeights,
elemCtx, m_simulator.model(),
m_element_chunks,
enableThreadParallel);
// Copy CPU vector to GPU vector using main stream and asynchronous transfer
m_weights->copyFromHostAsync(m_cpuWeights);
return *m_weights;
};
} else if (weightsType == "trueimpesanalytic") {
// Create CPU vector for the weights and initialize GPU vector
m_cpuWeights.resize(m_matrix->N());
m_pinnedWeightsMemory = std::make_unique<PinnedMemoryHolder<real_type>>(
const_cast<real_type*>(&m_cpuWeights[0][0]), m_cpuWeights.dim());
m_weights.emplace(m_cpuWeights);
const bool enableThreadParallel = m_parameters.cpr_weights_thread_parallel_;
// CPU implementation wrapped for GPU
weightsCalculator = [this, enableThreadParallel]() -> GPUVector& {
// Use the CPU implementation to calculate the weights
ElementContext elemCtx(m_simulator);
Amg::getTrueImpesWeightsAnalytic(pressureIndex,
m_cpuWeights,
elemCtx, m_simulator.model(),
m_element_chunks,
enableThreadParallel);
// Copy CPU vector to GPU vector using main stream and asynchronous transfer
m_weights->copyFromHostAsync(m_cpuWeights);
return *m_weights;
};
} else {
OPM_THROW(std::invalid_argument,
"Weights type " + weightsType
+ " not implemented for cpr."
" Please use quasiimpes, trueimpes or trueimpesanalytic.");
}
}
return weightsCalculator;
}
void updateMatrix(const Matrix& M)
{
if (!m_matrix) {
m_matrix.reset(new auto(GPUMatrix::fromMatrix(M)));
m_pinnedMatrixMemory = std::make_unique<PinnedMemoryHolder<real_type>>(
const_cast<real_type*>(&M[0][0][0][0]),
M.nonzeroes() * M[0][0].N() * M[0][0].M()
);
std::function<GPUVector&()> weightsCalculator = getWeightsCalculator();
m_gpuSolver = std::make_unique<SolverType>(
*m_matrix, isParallel(), m_propertyTree, pressureIndex, weightsCalculator, m_forceSerial, m_comm.get());
} else {
m_matrix->updateNonzeroValues(M, true);
m_gpuSolver->update();
}
}
void updateRhs(const Vector& b)
{
if (!m_rhs) {
m_rhs = std::make_unique<GPUVector>(b);
m_pinnedRhsMemory = std::make_unique<PinnedMemoryHolder<real_type>>(
const_cast<real_type*>(&b[0][0]),
b.dim()
);
} else {
// copy from host to device using main stream and asynchronous transfer
m_rhs->copyFromHostAsync(b);
}
}
FlowLinearSolverParameters m_parameters;
const Simulator& m_simulator;
const bool m_forceSerial;
PropertyTree m_propertyTree;
ElementChunksType m_element_chunks;
int m_lastSeenIterations = 0;
int m_solveCount = 0;
std::unique_ptr<GPUMatrix> m_matrix;
std::unique_ptr<SolverType> m_gpuSolver;
std::unique_ptr<GPUVector> m_rhs;
std::unique_ptr<GPUVector> m_x;
std::unique_ptr<PinnedMemoryHolder<real_type>> m_pinnedMatrixMemory;
std::unique_ptr<PinnedMemoryHolder<real_type>> m_pinnedRhsMemory;
std::unique_ptr<PinnedMemoryHolder<real_type>> m_pinnedXMemory;
std::unique_ptr<PinnedMemoryHolder<real_type>> m_pinnedWeightsMemory;
Vector m_cpuWeights;
std::optional<GPUVector> m_weights;
std::optional<GPUVectorInt> m_diagonalIndices;
std::shared_ptr<CommunicationType> m_comm;
std::vector<int> m_interiorRows;
std::vector<int> m_overlapRows;
std::any m_parallelInformation;
};
} // namespace Opm::gpuistl
#endif
Opm::AbstractISTLSolver
Abstract interface for ISTL solvers.
Definition AbstractISTLSolver.hpp:45
Opm::AbstractISTLSolver::checkConvergence
static bool checkConvergence(const Dune::InverseOperatorResult &result, const FlowLinearSolverParameters &parameters)
Check the convergence of the linear solver.
Definition AbstractISTLSolver.hpp:194
Opm::PropertyTree::get
T get(const std::string &key) const
Retrieve property value given hierarchical property key.
Definition PropertyTree.cpp:59
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< real_type >::fromMatrix
static GpuSparseMatrixWrapper< real_type > fromMatrix(const MatrixType &matrix, bool copyNonZeroElementsDirectly=false)
Definition GpuSparseMatrixWrapper.hpp:179
Opm::gpuistl::GpuVector
Definition gpu_type_detection.hpp:30
Opm::gpuistl::ISTLSolverGPUISTL::isParallel
bool isParallel() const
Check if we are running in parallel mode.
Definition ISTLSolverGPUISTL.hpp:315
Opm::gpuistl::ISTLSolverGPUISTL::ISTLSolverGPUISTL
ISTLSolverGPUISTL(const Simulator &simulator, const FlowLinearSolverParameters &parameters, bool forceSerial=false)
Construct a system solver.
Definition ISTLSolverGPUISTL.hpp:97
Opm::gpuistl::ISTLSolverGPUISTL::setMatrix
void setMatrix(const SparseMatrixAdapter &) override
Set the matrix for the solver.
Definition ISTLSolverGPUISTL.hpp:241
Opm::gpuistl::ISTLSolverGPUISTL::prepare
void prepare(const SparseMatrixAdapter &M, Vector &b) override
Prepare the solver with the given matrix and right-hand side vector.
Definition ISTLSolverGPUISTL.hpp:183
Opm::gpuistl::ISTLSolverGPUISTL::comm
const CommunicationType * comm() const override
Get the communication object used by the solver.
Definition ISTLSolverGPUISTL.hpp:305
Opm::gpuistl::ISTLSolverGPUISTL::getSolveCount
int getSolveCount() const override
Get the count of how many times the solver has been called.
Definition ISTLSolverGPUISTL.hpp:327
Opm::gpuistl::ISTLSolverGPUISTL::setResidual
void setResidual(Vector &) override
Set the residual vector.
Definition ISTLSolverGPUISTL.hpp:215
Opm::gpuistl::ISTLSolverGPUISTL::getResidual
void getResidual(Vector &b) const override
Get the residual vector.
Definition ISTLSolverGPUISTL.hpp:228
Opm::gpuistl::ISTLSolverGPUISTL::solve
bool solve(Vector &x) override
Solve the system of linear equations Ax = b.
Definition ISTLSolverGPUISTL.hpp:257
Opm::gpuistl::ISTLSolverGPUISTL::ISTLSolverGPUISTL
ISTLSolverGPUISTL(const Simulator &simulator)
Construct a system solver.
Definition ISTLSolverGPUISTL.hpp:136
Opm::gpuistl::ISTLSolverGPUISTL::iterations
int iterations() const override
Definition ISTLSolverGPUISTL.hpp:297
Opm::gpuistl::ISTLSolverGPUISTL::eraseMatrix
void eraseMatrix() override
Signals that the memory for the matrix internally in the solver could be erased.
Definition ISTLSolverGPUISTL.hpp:146
Opm::gpuistl::ISTLSolverGPUISTL::numAvailableSolvers
int numAvailableSolvers() const override
Get the number of available solvers.
Definition ISTLSolverGPUISTL.hpp:169
Opm::gpuistl::ISTLSolverGPUISTL::setActiveSolver
void setActiveSolver(int num) override
Set the active solver by its index.
Definition ISTLSolverGPUISTL.hpp:157
Opm::gpuistl::ISTLSolverGPUISTL::prepare
void prepare(const Matrix &M, Vector &b) override
Prepare the solver with the given matrix and right-hand side vector.
Definition ISTLSolverGPUISTL.hpp:197
Opm::gpuistl::detail::FlexibleSolverWrapper
FlexibleSolverWrapper is compilational trick to reduce compile time overhead.
Definition FlexibleSolverWrapper.hpp:45
Opm::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:233
Opm::setupPropertyTree
PropertyTree setupPropertyTree(FlowLinearSolverParameters p, bool linearSolverMaxIterSet, bool linearSolverReductionSet)
Set up a property tree intended for FlexibleSolver by either reading the tree from a JSON file or cre...
Definition setupPropertyTree.cpp:182
Opm::Parameters::IsSet
bool IsSet(bool errorIfNotRegistered=true)
Returns true if a parameter has been specified at runtime, false otherwise.
Definition parametersystem.hpp:270
Opm::Parameters::Get
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition parametersystem.hpp:187
Opm::FlowLinearSolverParameters
This class carries all parameters for the NewtonIterationBlackoilInterleaved class.
Definition FlowLinearSolverParameters.hpp:98