Loading...
Searching...
No Matches
AmgxInterface.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_AMGX_INTERFACE_HPP
#define OPM_AMGX_INTERFACE_HPP
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
#include <opm/simulators/linalg/gpuistl/detail/gpu_safe_call.hpp>
#include <opm/simulators/linalg/gpuistl/detail/gpu_type_detection.hpp>
#include <amgx_c.h>
#include <cuda_runtime.h>
#include <fmt/core.h>
#include <string>
#include <vector>
namespace Opm::gpuistl
{
class AmgxError : public std::runtime_error
{
public:
explicit AmgxError(const std::string& msg)
: std::runtime_error(msg)
{
}
};
inline std::string
getAmgxErrorMessage(
AMGX_RC err, const std::string& expression, const std::string& file, const std::string& function, int line)
{
char amgx_err_msg[4096];
AMGX_get_error_string(err, amgx_err_msg, sizeof(amgx_err_msg));
return fmt::format(
"AMGX error in expression: {}\nError code: {}\nError message: {}\nLocation: {}:{} in function {}",
expression,
static_cast<int>(err),
amgx_err_msg,
file,
line,
function);
}
inline void
amgxSafeCall(AMGX_RC rc, const std::string& expression, const std::string& file, const std::string& function, int line)
{
if (rc != AMGX_RC_OK) {
throw AmgxError(getAmgxErrorMessage(rc, expression, file, function, line));
}
}
#define OPM_AMGX_SAFE_CALL(expr) ::Opm::gpuistl::amgxSafeCall((expr), #expr, __FILE__, __func__, __LINE__)
class AmgxInterface
{
public:
static void initialize()
{
OPM_AMGX_SAFE_CALL(AMGX_initialize());
}
static void finalize()
{
OPM_AMGX_SAFE_CALL(AMGX_finalize());
}
static AMGX_config_handle createConfig(const std::string& config_string)
{
AMGX_config_handle config;
OPM_AMGX_SAFE_CALL(AMGX_config_create(&config, config_string.c_str()));
return config;
}
static AMGX_resources_handle createResources(AMGX_config_handle config)
{
AMGX_resources_handle resources;
OPM_AMGX_SAFE_CALL(AMGX_resources_create_simple(&resources, config));
return resources;
}
static AMGX_solver_handle createSolver(AMGX_resources_handle resources, AMGX_Mode mode, AMGX_config_handle config)
{
AMGX_solver_handle solver;
OPM_AMGX_SAFE_CALL(AMGX_solver_create(&solver, resources, mode, config));
return solver;
}
static AMGX_matrix_handle createMatrix(AMGX_resources_handle resources, AMGX_Mode mode)
{
AMGX_matrix_handle matrix;
OPM_AMGX_SAFE_CALL(AMGX_matrix_create(&matrix, resources, mode));
return matrix;
}
static AMGX_vector_handle createVector(AMGX_resources_handle resources, AMGX_Mode mode)
{
AMGX_vector_handle vector;
OPM_AMGX_SAFE_CALL(AMGX_vector_create(&vector, resources, mode));
return vector;
}
static void destroyConfig(AMGX_config_handle config)
{
if (config) {
OPM_AMGX_SAFE_CALL(AMGX_config_destroy(config));
}
}
static void destroyResources(AMGX_resources_handle resources)
{
if (resources) {
OPM_AMGX_SAFE_CALL(AMGX_resources_destroy(resources));
}
}
static void destroySolver(AMGX_solver_handle solver)
{
if (solver) {
OPM_AMGX_SAFE_CALL(AMGX_solver_destroy(solver));
}
}
template <typename MatrixType>
static void destroyMatrix(AMGX_matrix_handle amgx_matrix, const MatrixType& matrix)
{
if constexpr (!is_gpu_type<MatrixType>::value) {
// Unpin memory for CPU matrices
using T = typename MatrixType::field_type;
const T* values = &(matrix[0][0][0][0]);
OPM_AMGX_SAFE_CALL(AMGX_unpin_memory(const_cast<T*>(values)));
}
if (amgx_matrix) {
OPM_AMGX_SAFE_CALL(AMGX_matrix_destroy(amgx_matrix));
}
}
static void destroyVector(AMGX_vector_handle vector)
{
if (vector) {
OPM_AMGX_SAFE_CALL(AMGX_vector_destroy(vector));
}
}
template <typename T>
static void updateAmgxFromGpuVector(const GpuVector<T>& gpu_vec, AMGX_vector_handle amgx_vec)
{
// Get vector size from AMGX to verify compatibility
int n, block_dim;
OPM_AMGX_SAFE_CALL(AMGX_vector_get_size(amgx_vec, &n, &block_dim));
if (n > 0 && static_cast<size_t>(n * block_dim) != gpu_vec.dim()) {
throw AmgxError(fmt::format("Vector size mismatch in updateAmgxFromGpuVector: "
"AMGX vector size {} vs. GpuVector size {}",
n * block_dim,
gpu_vec.dim()));
}
// Get raw device pointer directly from GpuVector
const T* device_ptr = gpu_vec.data();
// Update AMGX vector directly with the device pointer
OPM_AMGX_SAFE_CALL(AMGX_vector_upload(amgx_vec, n, block_dim, const_cast<T*>(device_ptr)));
}
template <typename T>
static void updateGpuVectorFromAmgx(AMGX_vector_handle amgx_vec, GpuVector<T>& gpu_vec)
{
// Get vector size from AMGX to verify compatibility
int n, block_dim;
OPM_AMGX_SAFE_CALL(AMGX_vector_get_size(amgx_vec, &n, &block_dim));
if (static_cast<size_t>(n * block_dim) != gpu_vec.dim()) {
throw AmgxError(fmt::format("Vector size mismatch in updateGpuVectorFromAmgx: "
"AMGX vector size {} vs. GpuVector size {}",
n * block_dim,
gpu_vec.dim()));
}
// Get destination device pointer from GpuVector
T* dst_device_ptr = gpu_vec.data();
// Download data directly from AMGX vector to GpuVector's device memory
OPM_AMGX_SAFE_CALL(AMGX_vector_download(amgx_vec, dst_device_ptr));
}
template <typename VectorType>
static void transferVectorToAmgx(const VectorType& vec, AMGX_vector_handle amgx_vec)
{
if constexpr (is_gpu_type<VectorType>::value) {
// GPU implementation - device-to-device transfer
updateAmgxFromGpuVector(vec, amgx_vec);
} else {
// CPU implementation - host-to-device transfer
// For BlockVector, get dimensions
const int N = vec.size();
const int block_size = 1; // Assumes scalar vectors for CPU case
// Upload vector to AMGX
OPM_AMGX_SAFE_CALL(AMGX_vector_upload(amgx_vec, N, block_size, &vec[0][0]));
}
}
template <typename VectorType>
static void transferVectorFromAmgx(AMGX_vector_handle amgx_vec, VectorType& vec)
{
if constexpr (is_gpu_type<VectorType>::value) {
// GPU implementation - device-to-device transfer
updateGpuVectorFromAmgx(amgx_vec, vec);
} else {
// CPU implementation - device-to-host transfer
OPM_AMGX_SAFE_CALL(AMGX_vector_download(amgx_vec, &vec[0][0]));
}
}
template <typename T>
static void updateAmgxMatrixFromGpuSparseMatrix(const GpuSparseMatrixWrapper<T>& gpuSparseMatrix,
AMGX_matrix_handle amgxMatrix)
{
// Get matrix dimensions and sparsity information
auto n = detail::to_int(gpuSparseMatrix.N());
auto nnz = detail::to_int(gpuSparseMatrix.nonzeroes());
auto block_size = detail::to_int(gpuSparseMatrix.blockSize());
// Get device pointers directly from GpuSparseMatrix
const T* values = gpuSparseMatrix.getNonZeroValues().data();
const int* row_ptrs = gpuSparseMatrix.getRowIndices().data();
const int* col_indices = gpuSparseMatrix.getColumnIndices().data();
// Update AMGX matrix with the device pointers
OPM_AMGX_SAFE_CALL(
AMGX_matrix_upload_all(amgxMatrix, n, nnz, block_size, block_size, row_ptrs, col_indices, values, nullptr));
}
template <typename T>
static void updateAmgxMatrixCoefficientsFromGpuSparseMatrix(const GpuSparseMatrixWrapper<T>& gpuSparseMatrix,
AMGX_matrix_handle amgxMatrix)
{
// Get matrix dimensions and sparsity information
auto n = detail::to_int(gpuSparseMatrix.N());
auto nnz = detail::to_int(gpuSparseMatrix.nonzeroes());
// Get device pointer to values directly from GpuSparseMatrix
const T* values = gpuSparseMatrix.getNonZeroValues().data();
// Update coefficients in AMGX matrix with the device pointer
OPM_AMGX_SAFE_CALL(AMGX_matrix_replace_coefficients(amgxMatrix, n, nnz, values, nullptr));
}
template <typename T>
static void updateGpuSparseMatrixFromAmgxMatrix(AMGX_matrix_handle amgxMatrix, GpuSparseMatrixWrapper<T>& gpuSparseMatrix)
{
// Get matrix dimensions from AMGX
int n, nnz, block_sizex, block_sizey;
OPM_AMGX_SAFE_CALL(AMGX_matrix_get_size(amgxMatrix, &n, &block_sizex, &block_sizey));
OPM_AMGX_SAFE_CALL(AMGX_matrix_get_nnz(amgxMatrix, &nnz));
// Allocate temporary device memory for row pointers and column indices (required by AMGX API)
int* temp_row_ptrs;
int* temp_col_indices;
OPM_GPU_SAFE_CALL(cudaMalloc(&temp_row_ptrs, (n + 1) * sizeof(int)));
OPM_GPU_SAFE_CALL(cudaMalloc(&temp_col_indices, nnz * sizeof(int)));
// Get device pointer to values in the GpuSparseMatrix
T* gpu_values = gpuSparseMatrix.getNonZeroValues().data();
// AMGX requires a valid pointer for diagonal data (even if unused)
void* diag_data_ptr = nullptr;
// Download matrix values directly from AMGX to GpuSparseMatrix
OPM_AMGX_SAFE_CALL(AMGX_matrix_download_all(amgxMatrix,
temp_row_ptrs,
temp_col_indices,
gpu_values,
&diag_data_ptr));
// Clean up temporary device memory
OPM_GPU_SAFE_CALL(cudaFree(temp_row_ptrs));
OPM_GPU_SAFE_CALL(cudaFree(temp_col_indices));
}
template <typename MatrixType>
static void initializeMatrix(const MatrixType& matrix, AMGX_matrix_handle amgx_matrix)
{
if constexpr (is_gpu_type<MatrixType>::value) {
// GPU implementation - device-to-device transfer
updateAmgxMatrixFromGpuSparseMatrix(matrix, amgx_matrix);
} else {
// CPU implementation - host-to-device transfer
auto N = detail::to_int(matrix.N());
auto nnz = detail::to_int(matrix.nonzeroes());
const int block_size = 1; // Assumes scalar matrices for CPU case
// Setup sparsity pattern
std::vector<int> row_ptrs(N + 1);
std::vector<int> col_indices(nnz);
int pos = 0;
row_ptrs[0] = 0;
for (auto row = matrix.begin(); row != matrix.end(); ++row) {
for (auto col = row->begin(); col != row->end(); ++col) {
col_indices[pos++] = col.index();
}
row_ptrs[row.index() + 1] = pos;
}
// Get values pointer
using T = typename MatrixType::field_type;
const T* values = &(matrix[0][0][0][0]);
// Indexing explanation:
// matrix[0] - First row of the matrix
// [0] - First block in that row
// [0] - First row within the 1x1 block
// [0] - First column within the 1x1 block
// Pin memory for CPU-to-GPU transfers
OPM_AMGX_SAFE_CALL(AMGX_pin_memory(const_cast<T*>(values), sizeof(T) * nnz * block_size * block_size));
// Upload to AMGX
OPM_AMGX_SAFE_CALL(AMGX_matrix_upload_all(amgx_matrix,
N,
nnz,
block_size,
block_size,
row_ptrs.data(),
col_indices.data(),
const_cast<T*>(values),
nullptr));
}
}
static void initializeVector(int N, int block_size, AMGX_vector_handle amgx_vector)
{
OPM_AMGX_SAFE_CALL(AMGX_vector_set_zero(amgx_vector, N, block_size));
}
template <typename MatrixType>
static void updateMatrixValues(const MatrixType& matrix, AMGX_matrix_handle amgx_matrix)
{
if constexpr (is_gpu_type<MatrixType>::value) {
// GPU implementation - device-to-device transfer
updateAmgxMatrixCoefficientsFromGpuSparseMatrix(matrix, amgx_matrix);
} else {
// CPU implementation - host-to-device transfer
using T = typename MatrixType::field_type;
const T* values = &(matrix[0][0][0][0]);
OPM_AMGX_SAFE_CALL(AMGX_matrix_replace_coefficients(
amgx_matrix, matrix.N(), matrix.nonzeroes(), const_cast<T*>(values), nullptr));
}
}
template <typename MatrixFieldType, typename VectorFieldType>
static AMGX_Mode determineAmgxMode()
{
if constexpr (std::is_same_v<MatrixFieldType, double> && std::is_same_v<VectorFieldType, double>) {
return AMGX_mode_dDDI;
} else if constexpr (std::is_same_v<MatrixFieldType, float> && std::is_same_v<VectorFieldType, double>) {
return AMGX_mode_dDFI;
} else if constexpr (std::is_same_v<MatrixFieldType, float> && std::is_same_v<VectorFieldType, float>) {
return AMGX_mode_dFFI;
} else {
OPM_THROW(std::runtime_error,
"Unsupported combination of matrix and vector types for AMGX: "
+ std::string(typeid(MatrixFieldType).name()) + " and "
+ std::string(typeid(VectorFieldType).name()));
}
}
};
} // namespace Opm::gpuistl
#endif // OPM_AMGX_INTERFACE_HPP
Opm::gpuistl::AmgxError
Exception class for AMGX errors.
Definition AmgxInterface.hpp:44
Opm::gpuistl::AmgxInterface
Unified interface for AMGX operations with both CPU and GPU data structures.
Definition AmgxInterface.hpp:122
Opm::gpuistl::AmgxInterface::updateAmgxMatrixCoefficientsFromGpuSparseMatrix
static void updateAmgxMatrixCoefficientsFromGpuSparseMatrix(const GpuSparseMatrixWrapper< T > &gpuSparseMatrix, AMGX_matrix_handle amgxMatrix)
Update only the coefficient values of an AMGX matrix from a GpuSparseMatrix.
Definition AmgxInterface.hpp:446
Opm::gpuistl::AmgxInterface::determineAmgxMode
static AMGX_Mode determineAmgxMode()
Determine the appropriate AMGX mode based on matrix and vector field types.
Definition AmgxInterface.hpp:613
Opm::gpuistl::AmgxInterface::updateGpuVectorFromAmgx
static void updateGpuVectorFromAmgx(AMGX_vector_handle amgx_vec, GpuVector< T > &gpu_vec)
Update a GpuVector from an AMGX vector (device-to-device transfer).
Definition AmgxInterface.hpp:340
Opm::gpuistl::AmgxInterface::destroyResources
static void destroyResources(AMGX_resources_handle resources)
Destroy an AMGX resources handle.
Definition AmgxInterface.hpp:244
Opm::gpuistl::AmgxInterface::createResources
static AMGX_resources_handle createResources(AMGX_config_handle config)
Create AMGX resources from a config.
Definition AmgxInterface.hpp:167
Opm::gpuistl::AmgxInterface::initialize
static void initialize()
Initialize the AMGX library.
Definition AmgxInterface.hpp:130
Opm::gpuistl::AmgxInterface::destroyVector
static void destroyVector(AMGX_vector_handle vector)
Destroy an AMGX vector handle.
Definition AmgxInterface.hpp:292
Opm::gpuistl::AmgxInterface::createSolver
static AMGX_solver_handle createSolver(AMGX_resources_handle resources, AMGX_Mode mode, AMGX_config_handle config)
Create an AMGX solver.
Definition AmgxInterface.hpp:183
Opm::gpuistl::AmgxInterface::updateMatrixValues
static void updateMatrixValues(const MatrixType &matrix, AMGX_matrix_handle amgx_matrix)
Update matrix values in AMGX.
Definition AmgxInterface.hpp:590
Opm::gpuistl::AmgxInterface::updateAmgxFromGpuVector
static void updateAmgxFromGpuVector(const GpuVector< T > &gpu_vec, AMGX_vector_handle amgx_vec)
Update an AMGX vector from a GpuVector (device-to-device transfer).
Definition AmgxInterface.hpp:310
Opm::gpuistl::AmgxInterface::initializeVector
static void initializeVector(int N, int block_size, AMGX_vector_handle amgx_vector)
Initialize an AMGX vector with zeros.
Definition AmgxInterface.hpp:574
Opm::gpuistl::AmgxInterface::destroyConfig
static void destroyConfig(AMGX_config_handle config)
Destroy an AMGX config handle.
Definition AmgxInterface.hpp:226
Opm::gpuistl::AmgxInterface::createConfig
static AMGX_config_handle createConfig(const std::string &config_string)
Create an AMGX config handle from a configuration string.
Definition AmgxInterface.hpp:153
Opm::gpuistl::AmgxInterface::destroySolver
static void destroySolver(AMGX_solver_handle solver)
Destroy an AMGX solver handle.
Definition AmgxInterface.hpp:257
Opm::gpuistl::AmgxInterface::createMatrix
static AMGX_matrix_handle createMatrix(AMGX_resources_handle resources, AMGX_Mode mode)
Create an AMGX matrix.
Definition AmgxInterface.hpp:198
Opm::gpuistl::AmgxInterface::initializeMatrix
static void initializeMatrix(const MatrixType &matrix, AMGX_matrix_handle amgx_matrix)
Initialize an AMGX matrix from any matrix type (CPU or GPU).
Definition AmgxInterface.hpp:515
Opm::gpuistl::AmgxInterface::updateAmgxMatrixFromGpuSparseMatrix
static void updateAmgxMatrixFromGpuSparseMatrix(const GpuSparseMatrixWrapper< T > &gpuSparseMatrix, AMGX_matrix_handle amgxMatrix)
Update an AMGX matrix from a GpuSparseMatrixWrapper (device-to-device transfer).
Definition AmgxInterface.hpp:417
Opm::gpuistl::AmgxInterface::createVector
static AMGX_vector_handle createVector(AMGX_resources_handle resources, AMGX_Mode mode)
Create an AMGX vector.
Definition AmgxInterface.hpp:213
Opm::gpuistl::AmgxInterface::finalize
static void finalize()
Finalize the AMGX library.
Definition AmgxInterface.hpp:141
Opm::gpuistl::AmgxInterface::destroyMatrix
static void destroyMatrix(AMGX_matrix_handle amgx_matrix, const MatrixType &matrix)
Destroy an AMGX matrix handle.
Definition AmgxInterface.hpp:272
Opm::gpuistl::AmgxInterface::updateGpuSparseMatrixFromAmgxMatrix
static void updateGpuSparseMatrixFromAmgxMatrix(AMGX_matrix_handle amgxMatrix, GpuSparseMatrixWrapper< T > &gpuSparseMatrix)
Update a GpuSparseMatrixWrapper from an AMGX matrix (device-to-device transfer).
Definition AmgxInterface.hpp:474
Opm::gpuistl::AmgxInterface::transferVectorToAmgx
static void transferVectorToAmgx(const VectorType &vec, AMGX_vector_handle amgx_vec)
Transfer vector to AMGX from any vector type (CPU or GPU).
Definition AmgxInterface.hpp:370
Opm::gpuistl::AmgxInterface::transferVectorFromAmgx
static void transferVectorFromAmgx(AMGX_vector_handle amgx_vec, VectorType &vec)
Transfer vector from AMGX to any vector type (CPU or GPU).
Definition AmgxInterface.hpp:396
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::getRowIndices
GpuVector< int > & getRowIndices()
getRowIndices returns the row indices used to represent the BSR structure.
Definition GpuSparseMatrixWrapper.hpp:267
Opm::gpuistl::GpuSparseMatrixWrapper::getColumnIndices
GpuVector< int > & getColumnIndices()
getColumnIndices returns the column indices used to represent the BSR structure.
Definition GpuSparseMatrixWrapper.hpp:287
Opm::gpuistl::GpuSparseMatrixWrapper::getNonZeroValues
GpuVector< T > & getNonZeroValues()
getNonZeroValues returns the GPU vector containing the non-zero values (ordered by block)
Definition GpuSparseMatrixWrapper.hpp:247
Opm::gpuistl::GpuSparseMatrixWrapper::blockSize
size_t blockSize() const
blockSize size of the blocks
Definition GpuSparseMatrixWrapper.hpp:316
Opm::gpuistl::GpuSparseMatrixWrapper::nonzeroes
size_t nonzeroes() const
nonzeroes behaves as the Dune::BCRSMatrix::nonzeros() function and returns the number of non zero blo...
Definition GpuSparseMatrixWrapper.hpp:237
Opm::gpuistl::GpuSparseMatrixWrapper::N
size_t N() const
N returns the number of rows (which is equal to the number of columns).
Definition GpuSparseMatrixWrapper.hpp:228
Opm::gpuistl::GpuVector
Definition gpu_type_detection.hpp:30
Opm::gpuistl::detail::to_int
int to_int(std::size_t s)
to_int converts a (on most relevant platforms) 64 bits unsigned size_t to a signed 32 bits signed int
Definition safe_conversion.hpp:52
Opm::gpuistl::is_gpu_type
Type trait to detect if a type is a GPU type.
Definition gpu_type_detection.hpp:40